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SYNOPSIS 

This  report  describes  three  topics  considered 
during  the  second  quarter  of  Contract  AF19( 628) -5981. 

Three  three  topics  reported  on  are: 

a)  the  estimation  of  the  spectra  of  tran¬ 
sient  signals  with  additive  noise, 

b)  the  theoretical  aspects  of  clustering 
seismometers  in  a  large  aperture  seismic 
array, 

c)  the  results  of  tests  of  a  revised 
automatic  pP  phase  detection  technique. 

A  method  has  been  devised  to  estimate  the  energy 
density  spectra  of  transient  events  in  such  a  way  that 
the  variance  of  the  spectral  estimates  produced  by  addi¬ 
tive  noise  is  minimized.  . 

The  study  of  clusters  shows: 

1)  greater  spacing  between  sensors  would  im¬ 
prove  the  S/N  performance  of  the  arrays, 

2)  a  way  in  which  to  determine  the  amount  of 
clustering  possible  without  S/N  degrada¬ 
tion  depending  on  noise  conditions  at  the 
receiving  site. 
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Tests  with  the  revised  automatic  pP  detector 
have  shown  100$  correct  results  with  surface  focus 
events  (6),  83$  correct  results  for  earthquakes  (12), 
with  the  remaining  earthquakes  producing  anomalies. 


# 
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SECTION  I 

Introduction: 

This  is  the  Second  Quarterly  report  on  Contract 
AP19 (628) -5981 .  It  describes,  in  part,  the  studies  con¬ 
ducted  during  the  period  from  15  June  to  15  September. 

In  addition  to  the  special  task  mentioned  above, 
three  topics  were  considered  in  detail.  These  topics  are: 

a)  The  development  of  a  computation  technique 
for  the  estimation  of  the  spectra  of  transi¬ 
ents  corrupted  by  additive  noise  which  min¬ 
imizes  the  variance  in  the  spectral  estimate, 

b)  A  theoretical  examination  of  the  effects  of 
clustering  on  the  signal-to-noise  properties 
of  LASA's  for  a  variety  of  noise  conditions, 

c)  An  extension  and  expansion  of  the  earlier  work 
on  the  automatic  identification  of  the  pP 
phase  for  shallow  earthquakes. 

In  the  study  of  estimation  of  the  spectra  of 
transients  imbedded  in  noise,  a  "smoothing"  technique  was 
developed,  somewhat  analogous  to  that  of  Blackman  and  Tukey 
for  estimating  the  spectra  of  random  processes. 

The  cluster  study  results  indicate  that  larger 
clusters  than  those  currently  used  in  Montana  should  be 
used  in  order  to  increase  the  signal-to-noise  ratio  gain 
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of  an  array.  In  addition,  a  method  is  given  to  choose 
the  maximum  amount  of  clustering  which  can  be  used  with¬ 
out  degrading  performance.  This  limit  is  set  by  the  dis¬ 
tribution  of  local  and  teleseismic  noise  at  a  given  re¬ 
ceiving  site. 

The  revised  pP  detection  scheme  was  tested  with 
12  earthquakes  and  six  surface  events .  All  of  the  earth¬ 
quakes  studied  had  depth  picks  from  other  sources.  None 
of  the  surface  focus  events  indicated  a  pP  phase,  while 
ten  of  the  quakes  showed  depths  in  good  agreement  with  the 
previous  data.  One  of  the  remaining  quakes  showed  a  peak 
in  the  test  statistic  for  a  depth  appreciably  shallower 
than  the  reported  depth,  while  the  other  showed  three  peaks 
in  the  test  statistic,  one  of  which  corresponded  to  the 
reported  depth. 
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SECTION  II 


A.  Introduction 

In  the  last  quarterly  report,  it  was  noted  that 
microseismic  noise  was  independent  from  cluster  to  cluster 
in  the  Montana  MSA.  Further,  it  was  observed  that  signal 
correlation  was  quite  high  (mean  value  of  0.97)  across  the 
array.  As  a  result,  the  signal-to-noise  enhancement  which 
should  be  obtained  from  phased  summation  of  the  clusters  of 
the  array  is  very  nearly  equal  to  the  number  of  clusters 
summed.  Actual  experience  at  LAS A  bears  out  this  estimate, 
since  the  on-line  gain  achieved  is  approximately  13  dE, 
and  the  number  of  clusters  summed  is  21.  The  work  reported 
on  in  this  section  represents  an  extension  of  this  previous 
effort  to  include  the  effects  of  clustering  elements  of  an 
array  on  the  signal-to-noise  enhancement  which  can  be 
achieved. 

Specifically,  three  types  of  signal  processing  will 
be  discussed  in  conjunction  with  clusters  for  a  MSA. 

1.  Delayed  sum  (DS)  processing,  in  which  the 
element  outputs  are  added  with  each  output 
delayed  in  order  to  have  coherent  addition 
of  signals. 

2.  Weighted  delayed  sum  (WDS)  processing, 
which  is  the  same  as  DS  processing  with  the 
addition  that  each  element  output  is  given 
a  weight  that  is  inversely  proportional  to 
its  noise  power. 

3.  Filter  and  sum  (FS),  in  which  each  element 

is  filtered  before  summing  to  form  a  maximum  - 
likelihood  estimate  of  the  noise  fiell  in 
order  to  remove  as  much  noise  as  possible 
from  the  summed  output. 
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Further  details  of  these  processing  methods  are  pre¬ 
sented  elsewhere  (4). 

Results  of  experiments  with  these  three  processing 
techniques  obtained  by  UED  and  Lincoln  Labs.  (1-3)  will  be 
used  here  to  arrive  at  recommendations  for  signal  processing 
on  a  LASA.  Also,  some  theory  dealing  with  DS  processing  is 
developed.  Furthermore,  it  is  shown  that  a  seismic  array  is 
actually  a  diversity  system  so  that  much  of  the  current  tech¬ 
nical  literature  dealing  with  diversity  systems  may  be  applied 
to  seismic  arrays. 
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B.  A  Comparison  of  DS  Processing:  with  PS  Processing 

Experiments  performed  by  UED  on  actual  events  show 
that  for  the  existing  LASA,  PS  is  somewhat  better  than  DS 
(l).  Lincoln  Lab's  experimental  results  show  that  PS  for  the 
present  LASA  is  much  better  than  DS  (2).  This  is  true  even 
if  prefiltering  of  the  signals  is  used  to  remove  noise  that 
is  out  of  the  signal  band,  even  though  the  PS  processing  gain 
is  reduced.*  However,  present  computational  facilities  can¬ 
not  handle  on-line  FS  processing  of  LASA. 

The  poor  results  of  DS  processing  on  the  present 
LASA  are  to  be  expected  because  the  small  subarrays  maintain 
high  noise  correlation  between  the  elements,  preventing  power 
gains  anywhere  near  the  number  of  elements  in  the  array. 
Recently,  experiments  at  Lincoln  Labs  (3)  have  shown  that  by 
increasing  the  size  of  a  subarray  from  a  7  km  diameter  to  a 
22  km  diameter  will  raise  the  DS  gain  to  almost  the  same  level 
(about  3  dB  less)  as  the  PS  gain,  whereas  the  PS  gain  is  in¬ 
creased  only  slightly.  The  conclusion,  then,  is  that  FS  pro¬ 
cessing  is  not  worth  the  additional  complexity  and  time  that 
it  requires  when  it  performs  only  slightly  better  than  DS 
processing. 


*  The  overal]  gain  of  PS  processing  with  prefiltering  is 
actually  increased  over  that  of  PS  processing  alone,  but 
the  amount  of  gain  due  to  the  PS  processing  is  reduced. 
The  rest  of  the  gain  is  a  result  of  prefiltering. 
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C .  Seismic  Array  as  a  diversity  System 


In  this  section  it  is  shown  that  DS  and  '.IDS  process¬ 
ing  are  two  types  of  diversity  combining.  This  is  a  field 
that  has  recently  received  much  attention  in  the  technical 
literature,  particularly  with  respect  to  troposcatter  communi¬ 
cations  . 

By  the  term  diversity  system,  Brennan  (5)  refers 
to  "a  system  in  which  there  are  two  or  more  closely  similar 
copies  of  some  desired  signal."  With  m  such  copies 

f ^ ( t )  i=l , • . . , m  , 

each  copy  composed  of  signal  and  noise, 


fi(t)  =  +  ni  (t) 

a  general  linear  combination  will  be  considered 

m 

f(t)  =  ai  fi  (t-Tj. )  +  ...  +  anfm(t_Tm)  “ 

i=i 

Each  f^(t)  is  weighted  by  a  combining  coefficient  a^  and  de¬ 
layed  by  time  so  that  when  the  a^*s  and  t^'s  are  properly 
chosen  we  can  expect  f(t)  to  be  "better"  in  some  sense  than 
any  f^t) . 

There  are  three  common  linear  combination,  or  diver¬ 
sity  combining,  techniques.  In  all  three  the  t^'s  are  chosen 
so  that  the  signals,  the  s^’s  are  aligned  along  the  time  axis 
at  the  point  of  summation.  The  first  type  is  equal-gain  di¬ 
versity,  in  which  all  the  a^’s  are  equal.  This  is  the  same 
as  the  DS  processing  mentioned  previously. 

Secondly,  there  is  maximal-ratio  diversity.  In  this 
technique  the  weighting  factors,  the  a^'s,  are  chosen  to  max¬ 
imize  the  signal-to-noise  ratio  of  the  sum.  If  is  the 
signal  power  at  the  element,  and  P.~  the  noise  power  at 
the  ith  element,  then  the  optimum  weighting  factors  for  noise 
which  is  uncorrelated  between  elements  can  be  shown  to  be  (5) 
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The  signal-to-noise  ratio  of  the  sum,  then,  is  the  sum  of  the 
signal-to-noise  ratios.  This  type  of  diversity  "becomes  weighted 
delayed  sum  (WDS)  processing  for  the  case  where  all  signals, 
are  equal,  so  that  the  weighting  of  each  element  is  inversely 
proportional  to  its  noise  power.  The  equal-gain,  or  DS,  tech¬ 
nique  is  a  special  case  of  this  technique  when  all  element 
signal  powers  Pg^  are  equal  and  all  element  noise  powers  Pn^ 
are  equal. 

The  third  type  of  diversity  is  selector  diversity. 

In  this  case,  the  input  with  the  highest  signal-to-noise  ratio 
is  chosen,  and  all  the  others  are  removed.  In  other  words, 
if  element  ,j  has  the  highest  signal-to-noise  ratio,  the  weight¬ 
ing  factors  are: 


a.  «  °>  . 

'  i, 

The  problem  involved  with  using  maximal-ratio  diver¬ 
sity  or  selector  diversity  in  a  LASA  is  that  the  signal  ampli¬ 
tude,  as  well  as  the  noise  power,  must  be  known.  The  signal 
amplitude  at  any  element  is  not  constant  for  different  events, 
nor  is  it  constant  over  the  array  for  any  single  event,  (ex¬ 
periments  conducted  at  Lincoln  Labs  (3)  have  unearthed  events 
in  which  the  signal  amplitudes  differed  as  much  as  9  to  1 
over  the  LASA.)  The  difficulty  in  estimating  the  short  dura¬ 
tion  signal  accurately  at  each  element  with  a  noisy  background 
makes  the  use  of  these  two  diversity  techniques  too  difficult. 

Of  the  two  remaining  techniques,  DS  (equal-gain  com¬ 
bining)  and  WDS  (the  variation  on  maximal-ratio  combining), 

WDS  requires  that  an  estimate  of  the  noise  power  be  made  at 
each  element  and  used  there  as  its  weighting  factor.  This  means 
that  WDS  processing  is  more  complex  than  DS  processing.  Fur¬ 
ther  comparison  of  the  two  techniques  is  presented  in  Section  D . 
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P.  A  Comparison  of  PS  Processing  with  WPS  Processing 


Because  the  signal  amplitude  at  each  element  is  not 
measured  in  the  LASA  processing  one  must  use  a  processing 
technique  that  does  not  require  this  information .  In  its 
absence,  it  is  reasonable  to  assume  that  all  elements  have 
the  same  signal  amplitude.  Under  this  assumption,  and  the 
assumption  of  uncorrelated  noise  from  element  to  element, 

WPS  processing  is  the  optimum  combining  technique.  However, 
to  the  extent  that  the  equal  signal  assumption  does  not  hold 
for  an  event  WPS  processing  gains  begin  to  fall  off. 

Experiments  conducted  at  Lincoln  Labs  (3)  have  shown 
that  the  WPS  processing  gain  in  the  LASA  is  significantly 
better  than  the  PS  processing  gain  when  the  7  km.  subarrays 
are  used.  When  the  subarray  size  was  expanded  to  22  km., 
the  two  gains  both  increased  because  of  the  decrease  in  noise 
correlation  between  the  elements.  The  more  significant  gain 
increase  was  in  PS  processing,  so  that  for  the  larger  subarrays 
the  PS  processing  gain  was,  on  the  average,  within  1  dB  of  the 
WPS  processing  gain.  With  only  a  1  dB  difference,  WPS  does 
not  seem  to  be  worth  the  extra  complexity. 

For  the  case  of  a  LASA  with  large  subarrays  (large 
enough  so  that  the  noise  between  elements  can  be  assumed  to  be 
decorrelated) ,  we  can  derive  an  expression  for  the  loss  in 
gain  by  using  PS  processing  instead  of  WPS  processing.  For 
the  same  set  of  input  signals  and  noises  to  the  two  process¬ 
ing  schemes,  the  loss  L  is  simply  the  ratio  of  the  output 
signal-to-noise  ratio  of  PS  (SNR^g)  to  the  output  signal-to- 
noise  ratio  of  WPS  (SNR^o)  .  Assuming  the  signal  to  be  per¬ 
fectly  correlated  over  the  array,  we  have 

m 


( 


£ 

i-x 


2 


SNR 


PS 
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where  is  the  signal  power  at  the  ito  element 

is  the  noise  power  at  the  ith  element 
m  is  the  number  of  elements  in  the  array. 


Also 


fp  i 

?  tisi  >2 

SNRWI)S  -  ^  ’ 


( 


m  ■. 

E  r- 

i=  1  ^Ni 


Then  the  loss  is  given  by 

( 


9  111 

.2  s' 


L  = 


SNR 

SNR 


PS 

WPS 


m  ( — 

/,  Psi)'  . 
i=l  i=l  PNi 


Ni  >* 


We  can  consider  the  signal  power  across  the  array  as 
a  random  variable,  of  which  each  P  .  is  a  sample.  The  noise 

S  X 

power  can  be  regarded  in  the  same  manner,  with  its  samples  P^. 
With  E  |  |  indicating  expectation,  we  can  approximate 


^  i - 1  '  / - 1 

- mE 


m  j 

&  FPI 


mE  i  p~  i 

Ni 


^  PNi  “  mB  lPlfl  J 


ft  Vv 

i=l  ~;i 


/P  1 

mE  1  h 

^Ni  J 
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Assuming  Pg^  and  PH{  to  be  independent  random  variables,  we 


■Ni 


can  write 


E  |  =  E  Pg^  }  E  {  )  • 

*Ni  J  S1  J  ^Ni 


The  equation  for  the  loss  in  gain  then  becomes 

1 


L  » 


.rNii B  (4r  * 


a) 


or 


L  ^ 


m 


m 

<y.S.  P 


,?n 


iSL 


Note  that  the  above  approximation  is  exact  when  the 
signal  power  across  the  array  is  constant.  The  accuracy  of  the 
approximation  begins  to  deteriorate  as  the  samples  of  signal 
power  begin  to  spread.  The  approximation  itself,  however, 
depends  on  the  samples  of  noise  power  only. 

If  we  assign  a  reasonable  probability  density  to 
the  variable  Pjj.  ,  we  can  coi  put.<?  the  loss  L  of  equation  (1). 

Eor  simplicity  we  shall  use  a  symmetrical  triangular  density, 
shown  in  Figure  1  .  Computation  of  the  loss  1  gives 

2 

T  1  2h 


E  E  lPNiJ  (2n  +  w)  in  (!L-3)-  (2n-w)ln(—£' 

n-2' 

where  n  is  the  mean  of  the  density  and  w  is  its  width. 

This  loss  is  plotted  in  Figure  2  as  a  function  of  — ,  a 

parameter.  Note  that  even  at  =  .5,  the  lowest  reasonable  value 
allowed  by  the  assumed  distribution,  the  loss  is  only  1.4  dB. 
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e(pNi) 


2 

CD 


FIGURE  1  -  ASSUMED  PROBABILITY  DENSITY  OF  NOISE 
POWER  ACROSS  THE  ARRAY 
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^  _  Mean 
u  Spread 


FIGURE  2  -  LOSS  OF  DS  PROCESSING  BELOW  WDS  PROCESSING 

VS.  MEAN/SPREAD  PARAMETER  OF  DENSITY  IN  FIGURE  1 
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Within  the  accuracy  of  the  approximation  to  the  loss  in  gain, 
then,  WPS  processing  offers  negligible  improvement  over  DS 
processing  in  exchange  for  a  sizable  increase  in  processing 
complexity. 
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E.  Array  Gain  with  DS  Processing: 

The  primary  reason  for  building  an  array  is  to  im¬ 
prove  the  signal-to-noise  ratio  over  that  present  at  a  single 
element.  This  improvement,  or  gain,  should  he  written  in  a 
functional  form  to  determine  which  variables  affect  it.  In 
this  section  the  gain  equation  is  derived  based  on  certain  s 
stated  assumptions,  the  most  restrictive  being  that  the  signal 
amplitude  is  the  same  at  all  element  outputs  and  that  the 
noise  power  is  the  same  at  all  element  outputs.  Although  these 
assumptions  have  frequently  been  violated  in  actual  events 
recorded  at  the  LASA,  the  gain  derived  in  this  section  may  be 
used  as  a  first  approximation  to  what  is  expected  in  a  seismic 
array  and  as  a  basis  for  the  design  of  the  array  when  DS  pro¬ 
cessing  is  employed. 

In  an  array  that  uses  DS  processing  (all  element 

outputs  appropriately  delayed  and  then  summed),  the  array  gain 

may  be  easily  computed  as  a  function  of  the  number  of  elements 

in  the  array,  the  average  signal  cross-correlation  between  all 

pairs  of  elements,  and  the  average  noise  cross-correlation 

between  all  pairs  of  elements.  With  m  elements,  let  the  output 

of  the  ith\element  be  x.  ,  The  array  output  then  is  simply  the 

^  m 

sum  of  all  the  element  outputs, y  The  average  power  in  the 

array  output  is:  1=^ 
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If  we  assume  stationarity  of  the  outputs,  at  least  for  the 
duration  of  a  signal,  we  can  write 


_  1  rT 

xi  xj  =  te  27  J.T  *i(t)  xj  (t)  dt' 


where  R^A  (t)  is  the  cross-correlation  function  of  x.  and  x ^ 
Under  the  assumption  that  the  output  power  is  the  same  at 

each  element  (s^  =  x) ,  we  can  write 


~2 


Xi  X3  “  RijCo)  "  Pij  X  ' 

where  is  the  correlation  coefficient  between  the 

outputs  of  the  1th  and  elements.  Thus, 


in  m  ■  ■  7* 

P  =  V  T.  p.  .  xd  =  (in  -  m)"?  x  +  mx  , 
i=l  d=i 


ij 


(2) 


where  "p  is  the  average  correlation  coefficient  of  all  possible 
element  pairs . 

If  we  now  consider  th  element  outputs  to  be  composed 

of  two  uncorrelated  components,  signal  and  noise,  the  average 

power  of  the  array  output  is  the  sum  of  the  signal  power  P_ 

s 

and  noise  power  Pjj,  so  that  applying  equation  (2)  gives 

2  —  ~~2  ~~2 

P_  =  (m  -  m)  p_  s  +  ms 

S 

2  —  ~7l  "TT 

and  P^.  =  (m  -  m)  Pn  n  +  m  n  , 

~2 

where  s  is  the  signal  power  at  each  element 


n  is  the  noise  power  at  each  element 
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P_  is  the  average  signal  correlation  between  all 
s 

element  pairs 

Pjj  is  the  average  noise  correlation  between  all 
element  pairs 

The  gain  in  signal- to-noise  ratio,  then,  is 

rs/P„  (m-1)  7-  +1 

G  =  - JL  = - —2. -  .  (3) 

5?  /  ^  (m-1)  PH  +1 

The  above  equation  can  be  applied  to  any  summed 
array  where  the  signal  amplitude  at  each  element  is  constant, 
the  noise  power  at  each  element  is  constant,  and  the  signal  is 
uncorrelated  with  noise*  In  the  case  of  a  seismic  array  like 
the  LASA,  however,  it  is  useful  to  break  the  noise  into  its 
two  contributors s  local  noise  and  teleseismic  noise. 

These  two  components  differ,  insofar  as  the  array 
is  concerned,  by  the  fact  that  the  wavelengths  in  the  plane 
of  the  array  of  the  local  noise  are  much  smaller  than  those 
of  the  teleseismic  noise.  Thus,  one  would  expect  the  local 
noise  correlation  to  fall  off  more  quickly  with  distance  than 
that  of  the  teleseismic  noise.  A  second  difference  is  that, 
if  the  noise  power  is  normalized  to  a  signal  a  given  distance 
away,  the  normalized  local  noise  power  varies  more  from  site 
to  site  than  does  the  normalized  teleseismic  noise  power. 

Because  of  these  differences  between  local  noise  and 
teleseismic  noise,  a  gain  formula  will  be  developed  that  sep¬ 
arates  these  two  noise  components.  In  this  derivation  several 
assumptions  will  be  made* 

1.  Signal  amplitude  is  the  same  at  all  seismometer 
outputs . 

2.  Local  noise  power  is  the  same  at  all  seismometer 
outputs . 

3.  Teleseismic  noise  power  is  the  same  at  all  seis¬ 
mometer  outputs. 
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4.  Signal,  local  noise,  and  teleseismic  noise  are 
all  uncorrelated. 

5o  The  average  cross-correlation  for  local  noise 
over  any  cluster  is  constant  for  all  clusters, 
but  local  noise  is  uncorrelated  between  clusters. 

6.  The  average  cross-correlation  for  teleseismic 
noise  over  any  cluster  is  constant  for  all 
clusters . 

7.  Signal  is  perfectly  correlated  between  all 
seismometers  (6). 

The  first  five  assumptions  are  essentially  the  same  as  were 
made  in  deriving  formula  (3)® 

Some  additional  definitions  are  needed  before  de¬ 
riving  the  gain  formula  in  its  two-component  noise  form. 

These  are: 

“2 

n^  =  local  noise  power  at  each  element 

T 

n,p  =  teleseismic  noise  power  at  each  element 

2  2 
R  =  nL  /  nT 

P^  =  local  noise  power  at  the  array  output 

Pm  =  teleseismic  noise  power  at  the  output  of  each 
cluster® 

Pm  =  teleseismic  noise  power  at  the  output  of  the 
array 

Pt  =  average  correlation  of  local  noise  across  any 
cluster. 

P*Ti  =  average  correlation  of  teleseismic  noise  across 
any  cluster 

"p,_2  =  average  correlation  of  teleseismic  noise  between 
clusters . 

c  =  number  of  elements  per  cluster 

m  =  number  of  clusters  in  the  array, 
c 

If  filtering  is  used  following  the  DS,  the  linearity 
of  the  technique  allows  us  to  picture  each  element  output  as 
being  prefiltered  before  summing.  The  noise  powers  and  corre¬ 
lations  considered  here  are  those  present  after  this  prefiltering. 
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Tne  signal  power  after  tne  phased  sum  is 

„  2  ~~2 
P  =  m  3  t 

s 

using  tne  assumption  tnat  the  signal  is  perfectly  correlated 
over  the  array.  The  local  noise  power  out  of  the  array,  based 
on  equation  (2)  and  the  assumed  zero  correlation  between 
clusters,  is 

P^  =  c  [(c-l)  Pj^  +  lj  n^  •  “  « 

The  teleseismic  noise  power  out  of  any  cluster  is 


pTc  C 


[  (c-1) 


Pfji^  +  1  rifp 


and  the  teleseismic  noise  power  out  of  the  entire  array  is 


PT  -  (i) 


(~  1)  p-9  +  1 


T2 


Tc 


giving 


P  a 
rT  c 


(  C  “1 )  P  rp-j^  +  1 


(m—c)  Pqj2  ^  ® 

the  gain  in  signal-to-noise  ratio  can  then  be  written  as 

T 

Pq 

G  = 


PL  +  PT 


“2 

nT 


nT  +  nL 


(nT  +  nL)  m 


[(c-l)PL  +  1 


n^  +  1/ c  f(m-c)  PT2  +  c]  [(c-1)pt1  +  l]  n 


T 


Using  R,  the  ratio  of  local  noise  to  teleseismic  noise  at  the 
site,  we  have: 

+  R)  m 


G  = 


[(c-l)PL  +  l]  R  +  ~  [(m-c)PT2  +  cj  [(c-1)pt1  +  1 


U) 
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The  importance  of  this  equation  can  be  seen  if  it  is 
applied  to  the  present  LASA,  where  c  =  25  and  m  =  525.  In  this 
case  the  clusters  are  small  enough  to  consider  the  teleseismic 
noise  to  be  perfectly  correlated  over  a  cluster,  PT1  =  1.  We 
shall  also  assume  that  the  clusters  are  separated  enough  so 
that  teleseismic  noise  is  uncorrelated  between  clusters,  =  0. 

Thus ,  for  LASA 


GLASA 


525_(1  +  R) 

(24Pl  +  1)  R  +  25 


(5) 


Equation  (5)  is  plotted  in  Figure  5  for  several 
values  of  R.  The  dotted  line  indicates  the  gain  achievable 
when  the  noise  between  elements  is  completely  uncorrelated. 

The  curves  show  the  degradation  in  gain  as  "p^  increases  and/ 
or  as  R  decreases.  Tne  degradation  in  gain  is  primarily  due 
to  the  correlation  of  teleseismic  noise  across  a  cluster,  and 
it  is  increased  by  increasing  amounts  of  local  noise  correla¬ 
tion  across  a  cluster.  The  effect  of  local  noise  correlation 
is,  of  course,  diminished  by  a  smaller  ratio  of  local  noise  to 
teleseismic  noise,  but  then  the  effect  of  teleseismic  noise 
correlation  is  enhanced. 

Clustering  of  seismometers  in  the  LASA  was  originally 
proposed  (7)  to  ease  signal  transmission  problems  over  the  en¬ 
tire  20C  km  diameter  array.  The  size  of  a  cluster  (7km.)  was 
chosen  so  it  would  reduce  local  noise,  which  nas  an  average 
wavelength  of  about  4  km.  Teleseismic  noise,  with  an  average 
wavelength  of  about  16  km.  is  well  correlated  over  the  7  km. 
cluster,  so  that  the  action  of  a  cluster  is  indeed  to  reduce 
only  local  noiBe.  Expanding  a  cluster  sufficiently  allows  it 
to  strongly  reduce  teleseismic  noise  as  well  as  local  noise. 
(This  is  supported  by  an  UED  experiment  (1)  in  which  the 
power  gain  of  an  array  formed  by  DS  processing  the  outputs  of 
one  element  from  each  LASA  subarray  was  found  to  be  the 
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FIGURE  3  -  DS  GAIN  OF  LASA  VS.  AVERAGE  LOCAL  NOISE  CORRELATION 

FOR  SEVERAL  RATIOS  OF  LOCAL  NOISE  TO  TELESEISMIC  NOISE 
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number  of  elements  .in  tnat  arrays,  indicating  no  noise  correla¬ 
tion  between  the  element: *)  The  action  of  an  expanded  cluster 
cannot  be  viewed  solely  as  suppression  of  local  noise.  The 
cluster  concept  now  becomes  elements  grouped  together  only  for 
cabling  efficiency* 

Since  an  expanded  subarray  reduces  teleseismic  noise 
as  well  as  local  noise*  at  a  quiet  site  such  as  the  Montana 
LASA  where  teleseismic  noise  is  dominant,  one  would  expect  to 
see  a  very  large  gain  improvement  in  DS  processing  by  using 
the  expanded  cluster*  A  Lincoln  Labs  experiment  (3)  shows 
that  this  is  indeed  the  case  when  a  subarray  diameter  of  22 
km,  instead  of  7  km.  was  used*  The  ratio  of  local  noise  to 
teleseismic  noise  at  an  array  site,  in  fact,  has  been  shown 
above  to  be  a  strong  factor  in  array  performance  with  DS  pro¬ 
cessing  and  noise  correlation  between  the  elements  so  that 
it  must  be  considered  in  array  design  and  performance* 

Considering  a  LASA  with  expanded  clusters,  we  see 
that  at  a  quiet  site  (R<1)  such  as  the  one  in  Montana,  further 
decorrelation  of  local  noise  by  expansion  of  cluster  size  does 
not  improve  the  gain  much*  However,  with  a  cluster  diameter 
of  22  km*  or  greater  the  teleseismic  noise  correlation  across 
a  cluster  will  decrease  appreciably  from  1,  providing  a 
large  gain  improvement  over  the  array  with  smaller  clusters. 
With  the  present  concentration  of  clusters  near  the  center  of 
the  Montana  LASA,  however,  as  the  cluster  size  is  expanded 
the  outer  elements  of  each  cluster  will  come  close  to  each 
other,  causing  an  increase  in  the  teleseismic  noise  correla¬ 
tion  between  clusters  and  a  corresponding  decrease  in  gain. 

To  take  full  advantage  of  the  additional  array  gain  provided 
by  cluster  expansion,  then,  the  dense  array  center  must  be 
thinned  by  removing  several  clusters  from  that  region  and 
placing  them  in  the  less  dense  outlying  regions  of  the  array. 

In  other  words,  a  more  uniform  distribution  of  clusters  is 
desired  * 
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The  relatively  high  density  of  elements  in  the  center 
of  the  array  might,  in  fact,  have  a  more  detrimental  effect  on 
the  array  gain,  an  effect  that  is  not  immediately  obvious.  It 
is  conceivable  that  thinning  the  center  of  the  array  by  re¬ 
moving  elements,  but  not  placing  them  anywhere  else  in  the 
array,  could  actually  increase  the  array  gain.  In  other  words, 
including  these  elements  in  the  array  not  only  wastes  seismo¬ 
meters,  but  also  reduces  the  array  gain. 

As  a  simple  example  of  a  case  in  which  the  addition 
of  an  element  can  actually  reduce  the  array  gain,  consider  an 
array  composed  of  three  equally  weighted  elements  that  are 
separated  so  tnat  there  is  no  noise  correlation  between  them 
(p^  s=  o).  If  the  signal  is  perfectly  correlated  (ps  =  l), 
then  application  of  equation  (3)  gives  a  gain  of  G=3.  If  a 
fourth  element  is  added  with  perfect  signal  correlation  near 
enough  to  one  of  the  original  elements  so  that  its  noise 
correlation  with  that  element  is  0.9,  but  its  noise  correlation 
with  the  remaining  two  elements  is  zero  (See  Figure  4),  the 
average  noise  correlation  becomes  =  .15.  The  gain  of  equa¬ 
tion  (3)  then  becomes 


G 


4 

3  ( .15)  +  1 


=  2.76, 


less  than  the  gain  G=3  with  the  three  element  array. 

A  general  condition  for  the  deterioration  of  the 
array  gain  by  the  addition  of  an  element  can  easily  be  de¬ 
rived.  Assuming  that  the  signal  is  perfectly  correlated  over 
the  array,  equation  (2)  gives  the  m-element  array  gain  as 

G  =  — F- _ _ 

m  1  +  (m-1)  P 

m 

Adding  an  element  gives 

r  -  m+l 

m+1  '  1  +  "  w 
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FIGURE  l* 


CORRELATIONS 

p23  =  p13  =  Pl^  =  P2H 
0.9 

SKETCH  OF  AN  ARRAY  IN  WHICH  AN  EXTRA 
ELEMENT  DECREASES  THE  GAIN 


NOISE 

p12 
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The  additional  element  causes  a  drop  in  gain 


G. 


when 


m+l  <  1 , 


m 


which  results  in 


1-P 


P_  .  i  -  PTO  =  AP> 
m+l  m 


m 


(6) 


m 


Equation  (5)  is  the  condition  for  loss  of  gain  by 
addition  of  an  element  to  the  array.  An  equivalent  condition 
may  be  derived  by  noting  that 


m  o 


pm  m(m-l)^2  ?=i  P;l»3 


where  P,  is  the  noise  correlation  between  the  1th  element 
1 » 3 

and  the  element.  Also 


m+l  (m+l Jm 


m  3 

I  I  P 

3=2  1=1 


m 


i,d 


£  Pi 

i=l  Xt 


m+l 


Using  equation  (5),  we  have 


2  '  io  ,1  m 

rm+TJm  /;  i.3  1 

0=2  i=l  i=l  1 ’ 


m+l 


m  o 

2  7  y  n 

m(m-l;  0=2  i=l  pi,o 


1  - 


m  o 

I  t  p 

m(m-l)  o=2  i=l  Pi,0  , 

- 2 - 

m 


which  results  in 


m 


&  ^ 


m+l  +  (2m+l)  (m-1)  p 


m 


m+l 


2m 


(7) 
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Equation  (7)  is  a  condition  for  loss  of  gain  with  an 
additional  element  that  puts  a  bound  on  the  sum  of  the  noise 
correlations  between  the  added  element  and  the  existing  ones  in 
terms  of  the  number  of  elements  and  the  average  noise  correla¬ 
tion  in  the  original  array.  Por  an  array  with  many  elements 
we  can  approximate  equation  (7)  by 


m 


(7a) 


If  the  array  already  has  fairly  high  noise  correlations  between 
its  elements  (such  that  mPm  »£)m  the  condition  for  a  decrease 
in  gain  with  an  additional  element  is  simply  that  the  average 
noise  correlation  between  the  additional  aliment  and  the  other 
array  elements  be  greater  than  the  average  noise  correlation 
across  the  original  array. 


As  far  as  the  1ASA  is  concerned,  this  detrimental 


situation  might  very  conceivably  occur  if  an  element  were  added 
in  the  dense  center  of  the  array,  where  its  average  correlation 
with  the  present  elements  would  be  relatively  high.  Conversely, 
the  DS  array  gain  might  very  conceivably  be  increased  by  re¬ 
moving  an  element  from  the  array  center,  and  it  would  be  fur¬ 
ther  increased  by  relocating  that  element  in  the  sparsely 
populated  outskirts  of  the  array.  Applying  this  reasoning  to 
more  than  one  of  the  central  elements  leads  to  the  requirement 
that  the  array  elements  be  distributed  more  evenly. 


-25- 


GENERAL  ATRONICS  CORPORATION 


F.  Conclusions 

Several  conclusions  can  be  drawn  from  this  report. 
The  first  is  that  the  present  7km  diameter  cluster  size  is 
too  small,  and  expansion  of  the  cluster  size  will  increase 
the  signal-to-noise  ratio  gain  of  the  LASA.  With  expanded 
clusters  DS  processing  gains  are  close  to  those  of  FS  pro¬ 
cessing  and  WDS  processing  without  the  additional  complexity 
of  these  two  techniques.  DS  processing  with  expanded 
clusters,  then,  is  the  type  of  processing  to  use. 

Along  with  cluster  expansion  it  is  also  desirable 
to  spread  the  clusters  far  enough  apart  so  that  there  is 
very  low  noise  correlation  between  elements  of  adjacent 
clusters,  because  noise  correlation  decreases  the  array  gain. 
It  is  shown  that  removal  of  an  element  whose  noise  is  rela¬ 
tively  highly  correlated  with  the  noise  at  other  elements 
may  actually  decrease  the  array  gain,  so  that  low  noise 
correlation  is  indeed  desirable.  The  combination  of  cluster 
expansion  with  cluster  separation  in  effect  means  a  more 
uniform  arrangement  of  elements  than  is  present  in  the 
existing  LASA. 

A  gain  formula  for  DS  processing  is  derived  that 
considers  the  local  noise  and  the  teleseismic  noise  separately 
because  the  two  have  different  spatial  variations  in  their 
correlations  and  different  power.  The  formula  is  given  by 
equation  (4). 

Finally,  it  is  shown  that  a  seismic  array  with  DS 
or  WDS  processing  may  be  viewed  as  a  diversity  system  with 
a  commonly  used  diversity  combining  technique.  The  benefit 
of  this  approach  is  that  it  can  take  advantage  of  the  im¬ 
portant  developments  and  theory  that  are  appearing  in  the 
technical  literature. 
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SECTION  III 


A .  Introduction 

It  has  been  suggested  that  earthquakes  and  nuclear 
blasts  result  in  seismograms  that  differ  significantly  in 
their  spectral  energy  content.  If  this  is  the  case,  the 
analysis  of  the  spectral  distribution  of  energy  might  provide 
useful  statistics  for  discriminating  between  these  two  events. 
This  chapter  is  concerned  with  the  estimation  of  the  energy 
density  spectrum  of  a  transient  signal  in  the  case  where  the 
desired  signal  is  corrupted  by  additive,  stationary,  Gaussian 
noise.  Although  this  investigation  is  motivated  by  the  poten¬ 
tial  seismic  applications,  it  is  felt  that  the  results  may  be 
of  interest  for  other  applications  as  well. 

There  is  a  considerable  literature  associated  with 
spectral  estimation,  but  the  majority  of  the  papers  are  con¬ 
cerned  with  the  problem  of  estimating  the  power  density  spec¬ 
trum  of  a  stationary  random  process  on  the  basis  of  one  sample 
function  observed  during  an  interval  of  limited  duration. 
Although  these  papers  are  relevant,  they  consider  a  problem 
that  is  fundamentally  different  from  the  one  of  interest  here 
—  that  of  estimating  the  energy  density  spectrum  of  a  transient 
which  is  corrupted  by  additive  stationary  noise.  We  are  aware 
of  two  papers,  by  Mclvor  [8]  and  by  Larrowe  and  Crabtree  [9], 
that  have  attempted  to  apply  a  concept  of  a  time-varying 
spectrum  to  transient  seismic  data;  however,  we  do  not  feel 
that  these  papers  have  handled  the  problem  satisfactorily. 

These  two  papers  will  be  discussed  in  a  later  section. 

The  problem  under  consideration  can  be  formulated  as 
follows.  There  is  a  finite  energy  "signal",  x(t),  for  which 
the  Fourier  Transform, 
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X(f)  =  /  x  ( t )  e  J2lrft  dt,  (1) 

—  oo 

is  assumed  to  exist.  It  is  desired  to  calculate  the  energy 
density  spectrum  (EDS),  |X(f)|L  for  this  signal.  The  available 
data  consist  of  a  record  of  x(t)+n(t),  where  n(t)  is  stationary 
Gaussian  noise,  with  autocorrelation  function, 

oo 

Rn(  t)  =  n(t  )n(t+  t)~  =  /  Sn(f)e-7lft  df,  (2) 

—  oo 

where  S  (f)  is  the  power  density  spectrum  (PDS)  of  the  noise. 

It  is  assumed  that  the  approximate  onset  time  of  x(t)  is  known, 
but  that  its  duration  is  somewhat  uncertain.  The  problem  is 
then  to  devise  a  method  of  estimating  |X(f) |  .  In  doing  this 

a  major  problem  will  be  the  variability  in  the  estimates 
caused  by  the  stochastic  nature  of  the  noise,  and,  in  trying 
to  reduce  this  variability,  the  resulting  distortion  in  the 
"signal"  contribution  to  the  estimate  will  have  to  be  considered. 

The  "direct"  method  of  estimating  spectra,  which  essen¬ 
tially  consists  of  calculating  the  magnitude-squared  of  the 
Fourier  Transform  of  a  finite-time  sample  of  the  process,  will 
first  be  considered.  It  will  be  observed  that  this  method 
has  serious  statistical  deficiencies  in  that  the  variance  of 
the  noise  contribution  to  the  estimate  is  always  greater  than 
the  square  of  the  expected  value  of  the  noise  contribution. 
Furthermore,  this  will  be  the  case  even  if  the  observed  sample 
is  multiplied  by  a  "time  window"  before  transforming.  One 
approach  to  alleviating  this  variability  problem  is  the  indirect 
method  discussed  at  length  by  Blackman  and  Tukey  [10].  This 
method  was  proposed  by  them  as  the  best  means  of  estimating 
the  PDS  of  a  stationary  process  and  it  has  since  been  widely 
applied  to  that  problem.  However,  as  will  be  discussed  below, 
it  does  not  appear  to  be  a  useful  method  of  estimating  the  EDS 
of  a  transient  corrupted  by  additive  noise.  From  these  dis¬ 
cussions  it  will  be  concluded  that  a  direct  method,  despite 
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its  inherent  statistical  shortomcings ,  is  the  appropriate  method 
for  this  problem.  Suggestions  will  also  be  made  for  improving 
the  statistical  quality  of  the  method  as  further  experimental 
data  become  available.  A  potentially  valuable  method  of  re¬ 
ducing  variability  will  then  be  discussed.  This  scheme  is  based 
on  subtracting,  from  the  observed  x(t)+n(t),  an  estimate  of  the 
noise  that  is  based  only  on  the  noise  preceding  the  onset  of  the 
signal.  The  last  section  consists  of  a  critical  discussion 
of  the  papers  of  Mclvor  [8]  and  of  Larrow  and  Crabtree  [9]. 
Finally,  details  of  derivations  are  presented  in  Appendix  A. 

There  are  a  few  topics  that  are  omitted  from  this  dis¬ 
cussion  that  would  have  to  be  considered  before  implementing 
these  recommendations  with  actual  data.  For  convenience,  con¬ 
tinuous  representations  of  the  waveforms  are  employed  through¬ 
out  the  report.  As  calculations  will  presumably  be  done  on 
digital  computers,  it  is  necessary  to  consider  the  problem  of 
discrete  representations  of  these  waveforms.  This  should  not 
be  particularly  difficult,  but  there  are  some  issues  to  be 
resolved.  In  applying  these  calculations  there  will  be  several 
parameters  to  be  chosen  and  varied  experimentally.  No  de¬ 
tailed  recommendations  for  choosing  and  varying  these  parameters 
are  included  in  this  chapter.  Finally,  there  are  potential 
problems  of  approximation  that  should  be  considered.  The  noise 
prediction  scheme  involves  operations  closely  related  to 
"whitening".  Some  consideration  should  be  given  to  the  sensi¬ 
tivity  of  the  overall  results  to  the  approximations  and  com¬ 
promises  that  will  of  necessity  be  involved  in  the  computer 
implementation  of  the  predictor. 
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B .  The  Time  Window 

It  is  assumed  that  all  calculations  are  performed  on 

z(t)  =  a(t) [x(t)+n(t) ]  (3) 

where  a(t)  is  called  the  "time  window",  and  is  assumed  to  have 
a  Fourier  Transform 

oo 

A(f)  =  /  a(t)  e_0  2irft  dt  (4) 

—  oo 

Specifically,  |Z(f)|^  will  be  calculated,  where 

oo 

Z(f)  =  /  a(  t )  [x  ( t )  +n(  t )  ]e-J' 2ir  *  dt 

—  oo 

=  zx(f)  +  Zn(f)  (5) 

Perhaps  the  simplest  choice  of  a(t)  would  be  1  for  the  time 
interval  that  is  thought  to  contain  "most  of"  the  signal,  and 
0  outside.  However,  there  are  two  reasons  for  choosing  a 
smoother  off-on  transition  for  a(t).  In  the  application  under 
consideration,  the  "duration"  of  x(t)  is  not  known.  Presumably, 
a(t)  should  be  "turned  off"  when  x(t)  becomes  small  compared  to 
the  noise  level.  An  intuitively  appealing  way  to  incorporate 
the  uncertainty  about  an  appropriate  turn-off  time  is  to  turn 
off  a(t)  gradually.  In  this  way,  times  with  good  signal-to- 
noise  ratios  would  be  weighted  more  heavily  than  those  without. 
The  signal  contribution  to  Z(f),  Z  (f)  is  simply 

zx(f)  =  A ( f )  &  X(f)  (6) 

where  &  indicates  convolution.  To  reduce  the  distortion  in 
going  from  X(f)  to  Z  (f),  A(f)  should  be  a  narrow  pulse  at 
the  origin  of  the  frequency  domain.  If  a(t)  is  a  rectangular 
pulse,  A(f)  will  have  relatively  large  "side  lobes".  These 
can  be  reduced  by  a  more  "rounded"  choice  of  a(t). 

In  the  discussion  which  follows, the  "width"  of  A(f) 
and  the  size  of  its  sidelobes  will  be  important  parameters. 

To  illustrate  the  relation  between  the  shapes  of  A(f)  and  a(t) 
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two  examples  are  presented  in  Fig.  5*  The  first  example  involves 
a  rectangular  shape  for  a(t),  which  leads  to  relatively  large 
sidelobes  in  A(f).  The  second  example  is  the  Hanning  Window, 
which  has  been  recommended  by  Blackman  and  Tukey  as  a  "lag 
window".  (The  lag  window  will  be  defined  and  discussed  below.) 

The  Hanning  window  is  more  rounded  than  the  rectangular  one  and, 
hence,  has  relatively  small  sidelobes  and  a  somewhat  wider  cen¬ 
tral  peak.  It  will  be  convenient  in  much  of  what  follows  to 
assume 

A  ( f )  =  0  | f |  >  1/T  (7) 

The  two  examples  just  considered  indicate  the  degree 
of  this  approximation. 

C .  Noise-Alone  Case 
1 .  Direct  Method 

In  order  to  characterize  the  variability  problems  that 
result  from  the  noise  it  is  easiest  to  begin  by  considering  the 
case  of  noise  alone,  i.e.,  the  case  of  x(t)  =  0.  'In  Section  D' 
these  results  will  be  generalized  to  include  the  signal  contribu¬ 
tion  and  any  interaction  terms. 

I  Zn(  f )  I  2  =  /dct/dea(a)a(B)n(a)n(6)e“J2lTf(a"6)  (8) 

Taking  expected  values,  and  interchanging  the  order  of  integra¬ 
tion  and  expectation, 

E[  |  Zn(  f )  |  2  ]  =  |  Zn(  f )  |  2  =  /dot/dga(ot)a(e)Rn(a-g)e-  27rf(a-B) 

(9) 

Recognizing  that  the  result  of  the  g  integration  is  simply 
a (ct)  k  Rn(a)e  J TT^'a 

I  Zn(  f )  I  2  =  Jdaa(a)  /  du  eJ27TUa  A(u)Sn(u+f)  (10) 
Interchanging  the  order  of  integration,  and  integrating  over  a. 
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FIGURE  5  -  EXAMPLES  OF  TWO  TIME-WINDOWS  AND  THEIR  TRANSFORMS 
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|Zn(f)|2  =  /du  A*(u)A(u)Sn(u+f) 

=  / dv  Sn(v) | A(v-f) |2 

=  Sn(f)  8  | A( f ) |2  (11) 

p 

Thus,  the  expected  value  of  |Z  (f)|  is  a  "smoothed"  version 
of  the  true  PDS,  Sn(f),  where  the  "smoothing"  is  accomplished 
by  a  convolution  in  the  frequency  domain. 

Similar  calculations,  which  employ  the  well-known 
property  of  zero-mean  Gauss  variables, 

n( t )n( t )n( a)n( B )  =  n( t )n( t )  n( a)n( B)  +  n(t)n(a)  n ( t ) n ( B ) 

+  n(t)n( B)  n( a )n( t)  (12) 


are  presented  in  Appendix  A,  andi-yield 

cov[ |Zn(f1) |2, |Zn(f2) |2]  =  |/duSn(u)A*(u-f1)A(u+f2) \c 


which  specializes  to 


+  | /duS  (u)A*(u-f, )A(u-f„) |2 

1  2  (13) 


var[  |Zn(f )  |2]  =  I /duSn  ( u)  |  A(  u-f )  |  2  |  2  +  |  JduSn  (  u )  A(  u+f )  A*  (  u-f  )  \  'c 

=  (E[ |Zn(f) |2])2  +  |/duSn(u)A(u+f)A*(u-f) \‘ 

(1*0 


This  last  result  is  a  generalization  of  the  standard  "periodo- 
gram"  result  (see,  e.g.,  Davenport  and  Root,  Cll],  pp.  107- 
108),  which  may  be  obtained  from  the  above  by  choosing 

A ( t )  =  1//T  0  <  t  <  T  (15) 

=  0  elsewhere 


The  important  observation  from  Equation  (14)  is  that,  no 
matter  what  a(t)  is  chosen 

var[ | Zn( f ) |2]-  (E[ |Zn(f) |2])2  (l6) 

In  particular,  using  longer  and  longer  observation  times  will 
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never  reduce  the  variance  to  less  than  the  square  of  the  mean. 
Thus,  | Zn ( f ) |  cannot  provide  a  "consistent"  estimate  of  S  (f). 

Equation  (l6)  implies  that  the  direct  method,  without 
any  modifications  or  further  calculations,  would  not  be  a  satis¬ 
factory  method  for  estimating  the  PDS  of  a  stationary  random 
process.  There  are,  of  course,  ways  to  modify  the  method  that 
would  improve  its  statistical  features,  and  these  provide  one 
way  of  motivating  the  indirect  method.  Before  considering  these 
modifications,  a  few  properties  of  Equations  (11),  (13)  and  (1*0 
should  be  noted. 

For  the  present  discussion,  it  is  assumed  that  a(t) 
is  zero  outside  of  an  interval  of  length  T.  It  is  further 
assumed  that  within  this  interval  a(t)  has  been  chosen  in  such 
a  way  as  to  minimize  the  width  of  A(f).  As  illustrated  by 
example  in  the  previous  section,  this  could  result  in  an  A(f) 
having  a  total  width  of  roughly  2/T.  Assuming,  therefore,  that 

A ( f )  =0  | f |  >  1/T  (17) 

Equation  (13)  implies  that  |Zn(f^)|^  and  |Z  ( f 2 ) | 2  are  uncor 
related  for 

| | f ! | - l f 2 | [  >  2/T  (18) 

i.e.,  for  frequencies  separated  by  more  than  2/T.  And,  Equa¬ 
tion  (1*0  implies  that 

var[ |Zn(f) |2]  =  (E[ |Zn(f) |2])2  |f|  >  1/T  (19) 

Considering  only  the  samples  | ( 2 K/T )  | 2  for  K  -  0, 
these  results  may  be  summarized  by 

I Zn^~ T^ I  =  ^2n^^  ®  I  A ( f )  |  Jf=(2K/T)  (20) 

c°v[|zn(^)  |2,  |Zn(f^)|2]  =  0  for  |K  |  / 1  £  |  (21) 
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var[ |Zn(|^) |2]  =  { 


Zn(2K/T) 


2|Zn(0) 


K  /  0 
K  =  0 


(22) 


2 .  Direct  Method  with  Further  Smoothing 

In  considering  these  equations  an  important  question  is 
how  rapidly  S  (f)  varies  with  f.  If  it  varies  slowly  compared 
to  2/T,  then  the  convolution  of  Equation  (20)  will  not  distort 
the  power  density  spectrum,  i.e. , 

Sn(f)  fl  | A ( f ) |2  =  cSn(f)  (23) 

If,  furthermore,  it  varies  slowly  compared  to  frequency  inter¬ 
vals  of  several  times  2/T,  then,  the  samples  |Z  (2K/T)|‘“,  will 
have  the  same  expectation  for  several  successive  values  of  K. 

In  this  case,  a  local  average  of  these  statistics  would  be  an 
appropriate  estimate  of  this  expectation  and,  since  successive 
values  of  | ( 2K/T ) |  are  uncorrelated,  the  local  average  would 
be  better  statistically  than  the  individual  samples.  A  conve¬ 
nient  measure  of  the  variability  of  a  single  statistic  is  the 
"coefficient  of  variation",  which  is  defined  as  the  ratio  of 
the  standard  deviation  of  the  mean.  It  is  easy  to  show  that 
adding  together  M  uncorrelated  random  variables  with  identical 
means  and  variances  results  in  a  new  random  variable  with  a 
coefficient  of  variation  that  is  1//M  times  that  of  the  original 
variable . 

To  summarize,  for  K>0,  |Z  (2K/T)|*  has  a  coefficient 
of  variation  of  unity.  If  T  is  large  enough  that  S  (f)  is 
essentially  constant  over  frequency  regions  of  length  2M/2T, 

p 

then  the  local  averages  taken  over  M  of  the  samples  |Zn(2K/T)| 
will  yield  spectral  estimates  with  coefficients  of  variation 
equal  to  1//M.  The  calculation  just  described  —  the  "direct" 
method,  followed  by  further  "smoothing"  in  the  frequency  domain 
—  could  be  a  perfectly  satisfactory  way  of  estimating  the  PDS 
of  a  stationary  random  process,  provided  it  were  possible  to 
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choose  T  large  enough  to  achieve  a  low  coefficient  of  variation 
for  the  spectral  estimates.  This  method  is  essentially  equiva¬ 
lent  to  the  "indirect"  method  of  Blackman  and  Tukey,  as  will  be 
discussed  below.  A  more  careful,  and  thorough  discussion  may, 
of  course,  be  found  in  Blackman  and  Tukey. 

An  operation  similar  to  that  of  adding  M  adjacent  values 
of  |Zn(2K/T)|2  would  be  convolving  | ( f ) | 2  with  Q(f),  where  Q(f) 
is  a  pulse  of  width  2M/T  located  at  the  origin  of  the  frequency 
domain.  Because  this  operation  does  in  a  sense  combine  M  uncor¬ 
related  samples,  one  would  expect  it  to  result  in  coefficient 
of  variation  of  roughly  1//M.  Defining  the  result  of  this  fur¬ 
ther  smoothing  as  P(f)  yields 

P(f)  =  |Zn(f )  |2  0  Q(f)  (2*1) 

and,  therefore, 

PTfT  =  |Zn(f)  |2  a  Q (  f ) 

•  sn(f)  a  | a ( f ) | 2  a  q ( f ) 

=  c  Sn(f)  a  Q(f)  (25) 

where  the  last  approximate  equality  is  based  on  |A(f)|  being 
narrow  compared  to  Q(f),  which  will  be  the  case  if  Q(f)  is 
sufficiently  wide  compared  to  2/T  that  a  low  coefficient  of 
variability  results.  P(f)  is  sometimes  called  a  "smoothed 
periodogram" . 


3 .  Indirect  Method 


To  relate  this  method  to  the  indirect  method  of  Blackman 
and  Tukey,  it  is  easiest  to  make  a  few  specific  assumptions  and 
definitions . 

Let 


a(  t ) 


1//T 

0 


0  <  t  <  T 
elsewhere 


(26) 


and 
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(27) 


C 


Then  it  is  easy  to  show  that 


°n(f)  =  /0n(l)  e  J2nfT  dT  =  lzn^f)|2 


(28) 


Let  D^( t)  be  a  symmetric  function  of  t  that  is  zero  for 
I t |  >  Tm,  and 


Q±(f)  =  / D± (  t  )  e  J2ltfT  dt 


(29) 


Then 


Pi(f)=lZn(f)  |2aQi(f)  =  <fn(f)SQi(f)  =  jrdT  e  :2lTf  TDi(  O0i(  t)  (30) 


This  last  equation  states  that  the  direct  method  fol¬ 


lowed  by  smoothing  in  the  frequency  domain  is  equivalent  to 
calculating  an  estimate  of  the  noise  autocorrelation  function, 

0n( t)  [Equation  27],  multiplying  this  by  a  "lag  window",  D^(t), 
and  then  calculating  the  Fourier  Transform  of  the  result.  By 
assumption  D^(x)  is  zero  for  |i|  >  T^.  T„  is  called  the  "maximum 

lag".  If  D, (t)  is  chosen  appropriately,  Q^(f)  will  have  a  width 
of  approximately  2/T^,  and  the  coefficients  of  variability  of 
the  resulting  spectral  estimates  will  therefore  be  approximately 


The  label  "indirect  method"  is  applied  to  the  Blackman- 


Tukey  approach  because  it  begins  with  the  estimation  of  the 
autocorrelation  function.  Actually,  Blackman  and  Tukey  recom¬ 
mend 


<6o(  t)  T-j  t  |  jfQ 


n( t)n(t+ | t | ) dt  =  T_ j t |  0n(i)  (31) 


as  the  estimate,  rather  than  0  ( t) .  C  ( t)  has  some  appeal 

n  oo 

as  the  "natural"  way  to  use  all  of  the  available  data  to  esti¬ 
mate  the  autocorrelation  function,  and  it  is  an  unbiased  esti¬ 
mate  in  that 


E[Co(  0]  =  Rn(t) 


(32) 
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However,  there  is  some  dispute  as  to  which  is  preferable,  and 
Parzen,  for  example,  recommends  the  biased  estimator  0  (t). 
According  to  Parzen  C 1 ^ 3 , 0n (  t )  has  a  smaller  mean  square  error 
than  CsjqCOj  and  the  Fourier  Transform  of  0n(x)  will  always  be 
positive  [see  Equation  28],  which  is  not  the  case  for  Cq^t). 

In  practice,  the  lag  window  is  applied  in  two  steps. 

The  first,  and  usually  the  most  significant  in  its  effect,  is 
setting  the  estimated  autocorrelation  function  to  zero  for  all 
values  of  the  argument  larger  than  the  maximum  lag.  The  second 
is  the  multiplication  of  the  remaining  function  by  the  lag 
window.  In  the  discrete  formulation  appropriate  for  digital 
computers  it  is  often  easier  to  implement  the  second  step  by 
an  appropriate  convolution  in  the  frequency  domain  after  trans¬ 
forming  the  autocorrelation  function  (which  has  already  been 
truncated  at  the  maximum  lag).  The  problem  of  selecting  a  lag 
window  is  discussed  at  some  length,  with  several  illustrations, 
by  Blackman  and  Tukey.  There  is  no  straightforward  method  of 
determining  a  "best"  window,  but  it  does  not  seem  to  matter 
much  which  of  several  "good"  possibilities  are  used. 

In  using  these  calculations  to  estimate  the  PDS  of  a 
stationary  random  process,  it  is  first  necessary  to  make  some 
guesses  about  the  rate  of  change  of  S  (f)  with  frequency  in 
order  to  select  the  maximum  width  of  Q(f)  that  could  be  toler¬ 
ated  without  distortion.  It  is  also  necessary  to  specify  the 
desired  coefficient  of  variation  for  the  estimate;  for  the  sake 
of  illustration  assume  this  is  1/3-  Then  T  must  be  chosen  so 
that  about  nine  intervals  of  2/T  are  included  within  the  width 
of  Q(f).  Specifying  the  width  of  Q(f)  is  equivalent  to  specifying 
the  maximum  lag,  T^.  To  achieve  a  coefficient  of  variation  of 
l/3j  T  must  be  at  least  as  large  as  9TM. 

To  summarize,  the  direct  method  by  itself  leads  to  spec¬ 
tral  estimates  with  a  coefficient  of  variation  of  at  least  unity 
and  an  expected  value, 
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•Zn(f) I2  =  Sn(f)  S  lA(f) I2  (33) 

And,  this  will  be  the  case  no  matter  what  time  window,  a(t), 
is  chosen.  To  reduce  variability  it  is  necessary  to  perform 
additional  smoothing,  after  calculating  |Z  (f)|  (or,  equiva¬ 
lently,  in  the  indirect  method,  by  multiplying  the  autocorrelation 
function  by  a  lag  window).  This  additional  smoothing  leads  to  a 
spectral  estimate  P(f),  with  expected  value, 

P(f)  =  Sn(f)  fl  |  A  (  f )  |  2  g  Q  ( f )  (3*0 

Typically  Q(f)  is  much  broader  than  lACf)^,  so  that 

| A( f ) I  2  a  Q(f)  =  cQ( f )  (35) 

and  hence 

PTTT  =  cSn(f)  a  Q(f)  (36) 

Furthermore,  the  coefficient  of  variability  of  P(f)  will  be 
proportional  to  the  square  root  of  the  ratio  of  the  widths  of 

p 

|  A ( f )  |  and  Q(f).  This  ratio  may  be  expressed  as  T.^/T,  where 
Tv  is  the  "maximum  lag"  and  T  is  the  total  observation  interval. 

D.  Signal  Plus  Noise  Case 
1 .  Direct  Method 

The  observed  waveform  z(t)  is  now  assumed  to  include 
the  transient  "signal",  x(t). 

z(t)  =  a(t ) [x(t)+n(t ) ]  (37) 

As  shown  in  Appendix  A,  this  leads  to 

E [ | Z ( f ) |2]  =  |  Z  x ( f )  |2  +  | Z  n ( f )  |2  (38) 

var[ | Z ( f ) | 2 ]  =  var[ |Zn(f) |2]  +  2 | Zx ( f ) | 2 | ZR( f ) | 2 

+  2ReZ  2(f) /duS  (u)A*(u+f)A(u-f)  (39) 

A  11 
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And,  as  before,  for  frequencies  f,  and  f ,  separated  by  more  than 

2  p 

the  width  of  A(f),  |Z(f^)|  and  |Z(f  ,)|  are  uncorrelated  random 

variables.  (In  Equations  38  and  39,  Z  (f)  is  as  defined  earlier 

<A 

in  Equation  6,  Z  (f)  is  the  quantity  discussed  in  Section  C  ajnd 

n  ^ 

the  mean  and  variance  of  |Z  (f)|  are  given  by  Equations  11  and 
14.)  By  examining  Equations  (38)  and  (39),  it  may  be  observed 
that  the  expected  value  of  the  estimate  is  just  the  sum  of  the 
contributions  of  signal  and  noise,  but  that  the  variance  includes 
an  interaction  term  proportional  to  both  |Zv(f)|£‘  and  |Z  (f)| 
plus  another  interaction  term  that  will  be  zero  for  f  larger  in 
magnitude  than  half  the  width  of  A(f).  The  implications  of 
these  results  will  be  discussed  after  considering  the  indirect 
method . 

2 .  Indirect  Method 

To  dismiss  the  possibility  of  applying  the  indirect 
method,  with  the  customary  choice  of  parameters,  to  the  problem 
of  estimating  the  EDS  of  a  transient,  it  is  sufficient  to  con¬ 
sider  the  effect  of  the  indirect  calculation  on  the  signal  term 
alone.  In  Section  C  the  case  of  noise  alone  was  considered. 

The  total  spectral  estimate  includes  both  of  these  plus  an 
interaction  term.  The  signal  contribution,  obtained  by  setting 
n(t )  =  0 ,  is 

PY(f)  =  | Z(f) | 2fiQ(f)  =  |X(f)flA(f) | 28Q(f)  (40) 

yv  A. 

where  a(t)  is  a  rectangular  pulse  of  duration  T,  and  A(f)  there¬ 
fore  has  its  first  zero  crossings  at  ±1/T,  and  a  width  of  approx¬ 
imately  2/T.  The  customary  choice  of  parameters  for  the  indirect 
method  leads  to  Q(f)  having  a  width  of  approximately  ten  times 
that  of  A(f).  In  this  application,  that  would  involve  smoothing 
over  a  frequency  interval  of  20/T,  which  could  easily  be  com¬ 
parable  to  the  frequency  interval  covered  by  the  entire  spectrum 
of  X(f).  In  order  to  view  the  problem  in  the  time  domain,  it 
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should  be  recalled  that  T  is  chosen  to  correspond  roughly  to 
the  duration  of  x(t).  A  pulse  of  length  T  has  an  autocorrela¬ 
tion  function  extending  from  -T  to  +T.  Applying  the  indirect 
method  with  =  (1/1QJT  corresponds  to  truncating  this  autocor¬ 
relation  function  at  1/10  of  its  duration.  In  short,  using 
the  indirect  method  with  parameters  chosen  to  constrain  the 
standard  deviation  of  the  spectral  estimate  to  be  a  small  frac¬ 
tion  of  the  expected  value  of  the  noise  contribution  to  the 
spectral  estimate,  results  in  a  potentially  severe  distortion 
of  the  signal's  EDS. 

E .  Recommended  Method 

The  recommended  method  is  essentially  the  direct  one, 
except  that  it  could  be  followed  by  further  smoothing  of  the 
kind  accomplished  by  the  indirect  method.  However,  the  parti¬ 
cular  smoothing  that  would  be  appropriate  depends  on  the  de¬ 
tailed  statistics  of  this  process,  and  cannot  at  this  time  be 
specified.  To  understand  the  considerations  that  should  go 
into  this  smoothing,  it  is  useful  to  review  the  properties  of 


the  direct  method. 

As  before, 

z(t)  =  a(t) [x(t)+n(t) ]  (4l) 

Defining 

yK  =  |  Z ( 2K/T )  |2  (42) 

the  earlier  results  may  be  summarized  as  follows : 

y^  =  |Zx(2K/T) |2  +|Zn(2K/T) |2  (43) 

where 

|Zx(f) |2  =  | A ( f )  fi  X(f) |2  (44) 

| Z n ( f ) | 7  =  Sn(f)  S  | A( f ) | 2  (45) 

cov(yK,y^)  =  0  for  |K|  /  | i |  (46) 
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var  yK 


The  term 


[ I Zn( 2K/T) I  2  +  2|ZX(2K/T)  |2|Zn(2K/T)  |2](1  +  6Kj0) 

(47) 


1  K=0 

K  ,  0  0  K^O 


(48) 


results  from  the  integral  /duS^ ( u) A* ( u+f ) A ( u-f )  in  Equations 
(14)  and  (39).  This  integral  is  zero  for  |f|  >  1/T  and  equals 
|Zn(0) |2  for  f=0. 

Defining 

•hK  =  |Zn(2K/T)  |2  (49) 

XK  =  I  Zx(2K/T) I  2  (50) 

gives 

yK  =  XK  +  nK 

var  yK  =  +2xKnK]  ( 1+6K  Q)  (52) 

In  the  application  under  consideration  the  may  be 
assumed  known.  In  practice  these  could  be  obtained  by  per¬ 
forming  the  (indirect)  Blackman-Tukey  calculation  on  a  long 
section  of  noise  preceding  the  signal.  In  this  calculation 
the  lag  window  would  be  chosen  so  that  the  expected  value  of 
the  noise  contribution  to  both  calculations  would  be  the  same. 
This  implies 

D1(x)  =  a(  x)  8  a(-x)  (53) 

or  equivalently 

Q1(f)  =  | AC  D |2  (54) 

The  remaining  statistical  problem  is  then,  given  the 
random  variables  y^,  and  assuming  the  are  known,  how  should 
the  x^  be  estimated?  An  interesting  aspect  of  this  estimation 
problem  is  that  not  only  the  mean  of  y^ ,  but  also  its  variance, 
depend  on  xK* 
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If  nK«xK  little  difficulty  would  be  expected  in  esti¬ 
mating  x y  ,  but  in  the  case  where  =  Xj.  >  the  problem  becomes 
difficult.  There  are  two  general  ways  to  proceed.  The  first 
is  to  consider  each  y^  separately  and  develop  some  rule  for 
estimating  x  from  it.  An  obvious  rule,  but  certainly  not  the 
only  possibility,  would  be  to  use  the  unbiased  estimator  y ^  —  n ^  , 
which,  of  course,  has  the  same  variances  as  y^,  ( +2x^nK)  (  1  +  <5k  /  q 
This  approach  would  give  good  results  where  the  signal-to-noise 
ratio  (perhaps  defined  as  X^/t^)  is  large,  but  might  not  be 
satisfactory  where  it  is  small.  The  other  approach  would  be 
to  try  to  improve  the  statistics  having  high  coefficients  of 
variability  by  in  some  way  combining  adjacent  statistics.  Since 
the  y^  are  uncorrelated,  rules  for  combining  that  would  lead 
to  a  reduced  coefficient  of  variability  in  the  new  statistic 
could  be  specified.  Of  course,  the  new  statistic  would  be  a 
measure  of  the  energy  spectrum  in  a  broader  region  of  frequencies, 
and  it  probably  would  not  be  simply  the  total  energy  in  this 
region  of  frequency. 

If  it  is  desired  to  combine  adjacent  statistics,  there 
are  again  two  alternatives.  The  first,  and  simplest,  is  simply 
to  add  them.  If  they  had  comparable  means  and  variances,  this 
would  result  in  a  new  statistic  with  a  smaller  coefficient  of 
variation  (by  a  factor  of  1//2  if  pairs  are  combined).  In 
general,  however,  it  would  be  better  statistically  to  take  a 
weighted  sum,  where  the  weighting  would  depend  on  both  the 
mean  and  variance  of  each  statistic.*  The  disadvantages  of  this 

*  For  example  if  it  were  desired  to  choose  a  and  6  in  such  a 
way  as  to  minimize  the  coefficient  of  variation  of  2  =  ax+gy, 
where  x  and  y  are  uncorrelated  random  variables,  a  suitable 
choice  would  be 

a  =  (1/x) [u  /(u  +u  ) ] , 
y  A  j 

and  similarly  for  g,  where 
u  =  varx/x 

X 

This  results  in  u  =  (u  u  )/(u  +u  ) 

z  x  V  XV 
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weighting  are  twofold.  To  choose  the  weights  it  is  necessary  to 
know,  or  make  some  assumptions  about,  both  and  n  ,  and  the 

statistics  that  result  from  a  weighted  sum  would  correspond 
roughly  to  some  unequal  weighting  of  different  frequencies  in 
the  EDS. 

For  the  present  it  seems  reasonable  to  defer  development 
of  rules  for  combining  adjacent  statistics  until  some  further 
information  concerning  the  detailed  statistics  of  these  variables 
is  obtained.  Once  this  is  available,  it  may  turn  out  that  some 
compromise  between  simple  summation  and  an  optimally  weighted  sum 
seems  appropriate.  It  may  also  be  the  case  that  signal-to-noise 
ratios  for  some  of  the  y ..  are  typically  good  enough  that  no  com¬ 
bining  seems  necessary,  while  for  others  some  combining  always 
seems  necessary. 

F .  Improved  Method  Based  on  Subtracting  an  Estimate  of  the 

Noise  Waveform 

Because  the  "signal"  is  probably  confined  to  a  relatively 
short  time  interval  (e.g.,  2  or  ^  seconds)  and  because  the  noise 
appears  to  be  correlated  over  times  comparable  to  this,  it  should 
be  possible  to  reduce  the  variability  of  the  spectral  estimate 
by  subtracting  from  the  observed  signal  an  estimate  of  the  noise 
that  is  based  on  the  noise  preceding  the  observation  interval. 
Mathematically,  this  turns  out  to  be  very  similar  to  the  stan¬ 
dard  minimum-mean-square-error  filtering  problem  (see,  e.g.. 

Bode  and  Shannon,  [l3))>  but  there  are  some  important  differ¬ 
ences.  Specifically,  the  optimum  predictor  for  this  applica¬ 
tion  will  not  be  time  invariant,  as  it  is  in  the  standard 
smoothing  or  prediction  problem. 

For  this  discussion,  it  is  convenient  to  assume  that 
the  observation  interval  is  (0,T).  The  problem  is  to  devise 
an  estimate  n(t),  which  is  based  on  n(t)  for  t<0.  Before  an 
optimum  estimate  can  be  formulated  it  is  necessary  to  specify 
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some  criteria.  One  possibility  would  be  to  choose  n(t)  in  such 
a  ways  as  to  minimize 

E[(n(t)-n(t) )2] 

for  each  t.  A  more  appropriate  criterion  might  be  to  minimize 

p 

var[|Z(f)|  ]  where  z(t)  is  now  defined  by 

z(t)  =  [x ( t )+n( t ) -n( t ) ]  a(t)  (55) 

Unfortunately  this  latter  criterion  is  not  only  difficult  to 
apply  analytically,  but  it  cannot  be  applied  without  knowledge 

p 

of  |Z  (f)|  .  A  related,  but  simpler  criterion  would  be  to 

minimize  | Zn_^( f ) |  2,  the  expected  value  of  the  noise  contri¬ 
bution  to  the  spectral  estimate.  In  the  stationary  case  (with¬ 
out  n(t))  the  square  of  this  term  was  one  contribution  to  the 
variance  (Equation  1^).  The  same  result  holds  here  as  well 
(see  Appendix  A  and  details  below).  This  criterion  will  be 
discussed  after  considering  the  simpler  one  of  minimizing 
E[ ( n( t ) -n( t ) ) 2 ] .  It  will  turn  out  that  both  criteria  lead  to 
the  same  estimate. 

-  "  p 

1 .  Choosing  n(t)  to  Minimize  (n(t)-n(t)) 

In  addition  to  assuming  that  n(t)  is  a  zero-mean, 
stationary  Gaussian  process,  it  is  now  assumed  that  Sn(f) 
is  such  that  a  realizable  whitener,  with  a  realizable  inverse, 
exists.  [For  a  discussion  of  the  restrictions  on  Sn(f)  implied 
by  this  assumption,  see  Davenport  and  Root,  1959,  Chapter  11]. 
In  equation  form,  it  is  assumed  that  h-'*‘(t)  and  h(t)  exist, 
where  both  are  impulse  responses  of  realizable  systems  (i.e., 
h_'*'(t)  =  h(t)  =  0  for  t<0),  and 

h-1(t)  &  h(t)  =  uQ(t)  (56) 

h(t)  S  h(-t)  =  Rn(t)  (57) 

or,  equivalently, 

| H ( f ) |2  =  Sn(f)  (58) 
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Continuing  with  these  assumptions,  it  is  useful  to  consider 
the  following  hypothetical  system. 


h  1(t) 

w(t)  : 

h(t ) 

realizable 

realizabL 

whitener 

inverse 

FIGURE  6 

The  output  of  this  hypothetical  system  is  n(t)  because  of 
Equation  (56).  w(t)  is  white  noise,  with 

w(t)w(t+ t)  =  uq(t)  (59) 

Since  n(t)  can  be  recovered  immediately  from  w(t)  with  the 
realizable  filter  h(t),  nothing  is  lost  by  basing  the  estimate, 
n(t),  on  w(t)  for  t<0  rather  than  on  n(t)  for  t<0,  and  some 
mathematical  convenience  is  gained. 

By  inspection  of  Figure  6 


(a) 

w( a)  =  / 

n  (  B  )  h  1 ( a- 6 )  dB 

(60) 

—  00 

(t) 

n( t)  =  / 

w( a)h( t-a)  da 

(61) 

—  oo 


where  the  limits  in  parentheses  are  redundant  due  to  the 
realizability  of  h(t)  and  h-1(t). 

Rewriting  Equation  (6l) 

0  (t) 

n  ( t )  =  /  w(a)h(t-a)dct  +  /  w(a)h(t-a)da  (62) 

—  OO  0 

Because  w(t)  is  white,  and  the  two  integrals  involve  nonover¬ 
lapping  sections  of  w(t),  these  two  integrals  define  independent 
(zero  mean,  Gaussian)  random  processes.  The  first  integral 
represents  the  contribution  to  n(t)  of  the  white  noise  preceding 
t=0 .  The  second  integral  is  that  of  the  white  noise  following 
t=0.  To  minimize  ( n( t ) -n ( t ) ) c ,  the  first  integral  should  be 
chosen  as  n( t ) . 
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0 

n(t)  =  /  w(  a)h(t-ct)da 

—  oo 

0  (a) 

=  /  dah(t-a)  /  dgh  (a-3)n(3) 

—  OO  —  00 

0  0  -1 

=  /  d 3  /  doh( t-a)h  1(ot-g)n(g)  (63) 

(3) 


With  this  choice,  the  error,  n  (t),  will  be 

(t) 

n  (t)  =  n(t)-n(t)  =  /  w(a)  h(t-a)da  (64) 

0 

and,  as  already  observed,  ne(t)  and  n(t)  will  be  independent, 
zero-mean,  nonstationary  Gaussian  processes. 

The  autocorrelation  function  of  ne(t)  will  be  needed 
in  the  next  section.  In  Appendix  A  it  is  shown  that 


whereas 


b(t,i)  =  ng(t)ne(  t)  = 


min( t , 
(0) 


t) 

h(v)h( v+ 1 1—  t | ) dv 

(65) 


R  (t-i)  =  n ( t )n( t+T )  =  /  h(v)h(v+  t-T  | )dv  (66) 
n  (0) 

2 .  Choosing  n(t)  to  minimize  E[|Z  (f)|*~] 

As  indicated  earlier,  another  criterion  of  some  in- 
terest  would  be  choosing  n(t)  to  minimize  E  [  |  (  f )  | 1:1  ]  .  Because 

n(t)  is  assumed  Gaussian,  the  optimum  (in  the  sense  of  mini¬ 
mizing  | Z  ( f ) | ^  predictor  may  be  written  as  a  linear  operation 
on  n( t ) ,  t  <0 . 

0 

n  ( t )  =  /  g(t,a)n(a)  da  (67) 

—  OO 

i 

The  problem  is  to  choose  g(t,a)  to  minimize  |Z  (f)|  .  It  is 

not  at  first  obvious  whether  or  not  this  can  be  done  indepen¬ 
dently  for  all  f,  so  for  the  present,  f  is  regarded  as  a  single 
fixed  number.  The  details  of  the  derivation  appear  in  Appendix 
A,  and  only  an  outline  of  the  important  steps  appears  below. 
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Since , 


n  (t)  =  n(t)  -  /  g(t,o)n(a)  da 


(68) 


|Z  (f)|2  =  J/dtdi  a(  t  )a(  t)  e  1  2  T/nf  (t  )ne(x  ) 

(69) 

is  a  functional  of  g(t,a).  Considering 

g( t  j  a )  =  gQ ( t  j  a )  +  eg£ ( t  j  a )  (70) 


it  is  desired  to  choose  gg(t,a)  in  such  a  way  that  |Z  (f)| 
will  be  a  minimum  at  e  =  0,  for  any  choice  of  g  (t,a).  Choosing 
gg(t,a)  to  satisfy 
0 

/  dagQ(t  ,a)Rn(  a-6)  =  Rri(t-6)  for  6<0 


n 


(7D 


does  insure  that 


It  iVf)|2 


=  0 


e  =  0 


and 


(72) 


3  e 


^  Iz%(f>l2 


>  0 


e  =  0 


(73) 


(74) 


The  integral  equation,  (71) j  is  solved  by 

°  , 

gQ(t,a)  =  /  d6h  ( 6-a )  h(t-B) 

(a) 

Inspection  of  Equation  (63)  indicates  that  gQ(t,a), 

as  specified  by  Equation  (7*0,  leads  to  the  same  n(t)  as 

discussed  in  the  previous  section. 

With  this  choice  of  gQ  ( t ,  q_)  and  hence  of  n(t)  it  is 

possible  to  write  expressions  for  |Z  (f)|  and  var[  |  Z  (f)p], 

He  1  % 

but  unfortunately,  since  n  (t)  is  a  nonstationary  process,  these 

0 

equations  are  not  easy  to  simplify. 


|Zn(f) r  =  F ( f , -f ) 


(75) 
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var[|Zn(f) |2]  =  |P(f„-f)|2  +  |F(f,+f)|2 


(76) 


where 


P(f1,f2)  =  // dtdxa(t)a(  t)  / 


min( t ,  t  ) 


0 


dvh( v) 


h( v+ I t-x I )e 


-J2ir(f1t+f2T  ) 


(77) 


This  can  be  simplified  somewhat  to  yield 

OO  OO  Q 

F(f»-f)  =  2 /  da/  dg/  dva( a ) a( a+g )h( v) h( v+g ) cos2n f g 

(0)  0  (0) 


oo  o°  ct 

F( f if )  =  2 /  da/  dg /  dva( a) a ( a+g ) h( v) h( v+g ) e~ *  2 rf( 6  +  ‘  a ' 

(0)  0  (0) 


(78) 
+  2  a ) 

(79) 


In  both  cases,  the  corresponding  expressions  for  the  stationary 
case,  (|Z  (f)|  without  using  n(t))  may  be  obtained  by 
changing  the  upper  limit  on  the  v  integration  from  a  to  +<*>. 

In  the  stationary  case  these  expressions  can  be  simplified 
in  terms  of  the  Fourier  Transforms  of  a(t)  and  h(t).  Unfor¬ 
tunately,  comparable  simplification  does  not  appear  possible 
in  the  nonstationary  case  considered  here,  unless  specific 
assumptions  about  a(t)  and  h(t)  are  made,  and  thus  far  even 
this  has  not  helped  much. 

In  Appendix  A  expressions  are  given  for  F(f,-f), 
assuming  a(t)  and  h(t)  are  rectangular  pulses  of  arbitrary 
length.  In  the  case  where  they  have  the  same  length,  these 
expressions  reduce  to 

|Zn(f) |2  =  l/2|Zn(f) |2  (80) 

In  the  stationary  case,  discussed  earlier,  it  is  usually 

2 

possible  to  neglect  |F(f,+f)|  ,  the  second  term  in  the  expres- 

p 

sion  for  var[ | Zn(f ) |  ] .  It  is  not  yet  clear  whether  or  not 
this  is  also  true  of  the  nonstationary  case.  Even  for  the 
special  case  of  rectangular  a(t)  and  h(t),we  have  as  yet  been 
unable  to  obtain  a  simple  expression  for  |F(f,+f)|  . 
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Physically,  the  reduction  in  the  variance  of  the  spec¬ 
tral  estimate  would  be  expected  to  depend  on  the  ratio  of  the 
durations  of  a(t)  and  h(t),  with  the  variance  being  quite  small 
(i.e.,  n(t),  a  very  good  estimate)  when  a(t)  is  short  compared 
to  h(t)  and  relatively  large  when  the  reverse  is  true.  (If  h(t) 

=  0  for  t>B,  n(t)  =  0  for  t>B). 

In  summary,  it  appears  possible  to  reduce  the  noise  con¬ 
tribution  to  the  spectral  estimate  by  subtracting  an  estimate 
of  the  noise  waveform  that  is  based  on  the  noise  preceding  the 
signal.  The  analytical  predictions  of  the  resulting  variability 
in  the  estimate  have  yet  to  be  worked  out,  and  will  probably 
require  some  simplifying  assumptions  and  approximations.  The 
details  of  implementing  some  approximation  to  this  calculation 
on  a  digital  computer  also  have  to  be  worked  out.  In  general, 
the  output  of  the  whitener  depends  on  the  infinite  past,  and 
this  dependence  would  have  to  be  truncated  in  any  implementa¬ 
tion.  Actually  the  finite-discrete  version  of  the  problem  is 
essentially  the  same  as  the  linear  mean  square  regression  prob¬ 
lem  [see,  for  example,  Cramdr,  195^,  pp .  302-305]  and  should  be 
straightforward  to  program. 

G .  Some  Comments  on  the  Papers  of  Mclvor  and  by  Larrowe  and 

Crabtree 

Both  of  these  papers  attempt  to  apply  to  seismic 
waveforms  a  concept  of  a  time-varying  power  spectrum.  The  basic 
concept  can  be  described  in  the  notation  of  this  chapter  by 
allowing  the  location  of  the  time  window,  a(t),  to  vary.  In 
this  notation,  then, 

z  (t)  =  a( t-t  n ) y ( t )  (81) 

zQ  u 

and 

|Z  (f)|2  =  |/a(t-tn)y(t)e"J2irftdf  |  (82) 

u0  u 
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y(t)  is  used  here  for  the  seismic  waveform  to  avoid  specifying 
whether  it  is  a  transient  (as  x(t)  was)  or  a  continuing  process 
(as  n(t)  was).  |Ztg(f)|  would  be  labelled  by  words  such  as 
instantaneous  power  spectrum. 

With  this  formulation,  several  questions  immediately 
arise,  most  of  which  stem  from  the  nontrivial  question  of  what 
is  the  meaning  of  |  Z-t q ( f )  |  ?  Is  it  an  estimate  of  a  power  (or 
energy)  density  spectrum  or  is  it  an  interesting  functional 
of  y(t)  in  its  own  right.  Whereas  the  power  density  spectrum 
of  a  stationary  random  process  and  the  energy  density  spectrum 
of  a  transient  both  have  appealing  "physical"  interpretations,* 
it  is  not  at  all  clear  what  interpretation  is  to  be  given  to 
|  Z-t q  ( f )  |  »  unless  assumptions  are  made  about  how  slowly  it 

varies  with  tg.  In  our  opinion  there  is  a  burden  on  the  user  of 
a  "time-varying  spectrum"  to  demonstrate  its  physical  rele¬ 
vance,  its  utility,  or,  preferably,  both,  we  do  not  believe 
this  is  satisfactorily  done  by  either  of  the  papers  considered. 

Mclvor  completely  neglects  the  problem  of  variability, 
which  can,  of  course,  be  important.  If,  for  example,  y(t) 
is  a  stationary  Gaussian  process,  then,  no  matter  what  a(t) 
is  used 

var [ | Z .  ( f ) | 2]  -  (EC | Z.  (f)|2])2  (83) 

0  z0 

Furthermore,  he  is  fundamentally  in  error  when  he  asserts 
that  Blackman  and  Tukey's  lag  window  is  just  another  time 
window  and  claims  that  their  indirect  method  is  simply  a 
special  case  of  his  (essentially  direct)  method. 

Larrowe  and  Crabtree  briefly  mention,  and  then  ignore, 

the  problem  of  variability  by  making  assumptions  that  imply  it 

does  not  exist.  They  divide  the  observation  interval  into 

*The  average  power  (or  the  total  energy)  of  the  output  of  a 
narrow  band  filter  at  frequency  fg  is  proportional  to  the  value 
of  the  power  (or  energy)  density  spectrum  evaluated  at  f q . 
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several  subintervals,  and,  for  each,  they  apply  the  equation 
(their  Equation  36) 


0  —  0  ,  ; 

n  n-1  t  , 
n-1 


to  estimate  the  PDS.  In  this  equation,  f(t)  is  the  seismic 
waveform.  They  assume  (page  16)  "each  [subinterval]  long  enough 
to  permit  a  reasonably  accurate  computation  of  p(u)"  [using 
this  equation].  As  discussed  in  Section  C  of  this  chapter 
(and  elsewhere),  no  matter  how  long  the  interval,  if  f(t) 
is  a  stationary  Gaussian  process. 


var[p(w  )  ]  -  (pU  )  T2 


(85) 


Two  pages  later,  in  a  derivation  based  on  this  assumption, 

they  let  At  =  t  -t  ,  go  to  zero! 

J  n  n-1  & 


It  turns  out  that  the  result  of  this  derivation  is 


essentially  correct  and  is  of  some  interest.  This  result  is 
that  an  equivalence  exists  between  the  time  window  and  a 
suitably  defined  frequency  window.  This  equivalence  is  also 
demonstrated  by  Mclvor,  but  his  derivation  is  a  little  diffi¬ 
cult  to  follow  because  it  involves  several  normalizations 
that  are  needed  later  in  his  paper.  As  this  equivalence  is 
fairly  easy  to  demonstrate  in  the  notation  of  th'..;  u  • " . 

a  brief  derivation  follows. 


Zfc  (fQ)  =  /  a(t-tQ)y(t)e 


00 


-J2*f0t 


dt 


(86) 


Considering  a  filter  with  impulse  response 


hf  (t)  =  a(-t)  cos  2  uf Qt 


(87) 


and  output  w(t),  yields 
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w(t)  =  y ( t )  8  h  (t) 

*0 

=  /y(  T)a(  i-t)  cos2iTf0(t- t)  di 
=  Re[/y(  x)a(  T-t)e'- “■  TT^°t  e  '^7T"°T 
=  Re  eJ2irfotZt(f0) 

=  |Zt(f0)  |cos(2irf0t  +<Zt(fQ))  (88) 

Thus  the  envelope  of  the  output  of  the  filter  at  time  t(  is 
equal  to  | Z-t q C f*Q )  |  »  and  therefore  using  the  time  window 
a(t-tg)  located  at  tg,  to  estimate  the  PDS  at  f(  is  in  a  sense 
equivalent  to  using  the  frequency  window  ^-[A#(  f-f g )+A* ( f+fg)  ]  , 
located  at  f  g ,  and  examining  the  envelope  at  time  t  .  Of 
course,  in  practical  situations,  calculating  the  envelope 
may  be  difficult,  but  in  the  case  where  |Ztg(fg)|  varies  only 
slowly  with  t  ,,  there  should  be  no  problem.  In  principle, 
another  filter  h  ,  ( t )  with 

hp(t)  =  a(-t)  sin2irfgt  (89) 

could  be  used,  in  which  case  the  sum  of  the  squares  of  the 
outputs  of  the  two  filters  at  time  tg  would  be  proportional 
to  |Zt0(f0)|2. 
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SECTION  IV 

Automatic  Toot  for  pP  in  Shallow  Earthquakes 

The  use  of  the  reflected  phase  pP  to  determine  the 
depth  of  shallow  earthquakes  is  often  rendered  ineffective  by 
the  masking  of  pP  by  the  coda  of  P.  Por  earthquakes  above  the 
Moho,  pP  is  usually  buried  in  the  first  13  seconds  of  coda 
following  P  onset  where  it  is  often  intractably  unrecognizable 
even  by  the  skilled  observer.  Consequently,  a  test  has  been  de¬ 
vised  to  extract  pP  from  the  overlying  coda  by  computer. 

The  method  is  based  on  signal  enhancement  of  pP 
through  phased  sums  of  a  seismic  array  of  continental  dimensions 
(3000-7000  kilometers).  This  processing  scheme  is  similar  to 
one  previously  reported*,  however,  several  important  additions 
have  been  made  to  facilitate  comparison  of  snallow  quakes  and 
surface  events.  (As  previously  given,  the  test  was  unable  to 
distinguish  between  the  two  classes.) 

Assume  a  set  of  stations  i=l,...,n  and  epicentral 
distances  a^.  Determine  the  angle  of  incidence  i^  for  each  a^ 
(Richter,  1958,  pp.  664-6).  Assume  a  set  of  test  depths,  h. 

o 

in  the  range  10  to  41  kilometers.  (Increments  of  1.6  km  were 
used).  Calculate  a  set  of  time  delays  of  P-pP  by: 

=  (2^/c)  cos  ^ 

Por  each  test  depth  h.  align  the  seismograms  for  proper  relative 

J 

displacement  according  to  t,  ..  Porm  a  weighted  sum  of  the  seis- 

**■  J 

mo grams,  by  multiplying  each  by  a  factor 


*GAC  Report  1456-2026-9,  January  1966. 
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S 


i 


{  (1/T-l) 

(1/C2) 


T1  2  I* 

j  si 

~Tp  2 

ni(t)dt 


where  and  N^  are  the  mean  energies  in  the  P  and  noise  regions 
respectively.  Por  this  computation  a  P  interval  of  one  second 
was  chosen.  This  interval  was  centered  about  the  first  definite 
zero  crossing  of  P.  The  noise  interval  chosen  was  from  13  to  3 
seconds  before  P  onset.  The  weighting  factor  is  derived  in 
Appendix  B,  using  the  constraint  that  the  signal-to-noi3e  ratio 
of  the  phase-sum  seismogram  be  a  maximum  with  respect  to  all 
other  possible  weighting  factors. 

After  alignment  and  weighting,  sum  seismograms  s.(t) 

D 

are  formed  for  each  test  depth  h^ .  A  window  of  duration  t  is 
placed  about  the  location  t .=(2hi/c)  in  each  s.(t),  and  the 
average  energy  in  that  window  measured: 


i  p*^  a  +T/2  p 

Sj  “  4  JtLt/2  SRt^t 

cl 


The  corresponding  average  energy  C.  ofthe  P-coda  of  S.(t)  is 

J  J 

measured  outside  the  window  about  t  . ,  and  a  test  statistic  formed 

J 


G  . 

z 


■  (VV 


A  plot  of  this  statistic  v.s.  h.  should  show  Tweaks  at  certain 

J 

values  of  h^.  One  reason,  in  the  case  of  a  shallow  earthquake, 

is  the  presence  of  pP  in  the  t.  interval.  A  second  reason  for 

J 

peaks  abserved  with  surface  and  deep  events  is  the  purely  ran¬ 
dom  occurrence  of  more  energy  in  one  interval  than  in  the  others. 
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The  ambiguity  is  easily  removed  by  formation  of  a  new  test  sta¬ 
tistic,  s.,  which  is  the  product  of  e.  and  the  average  cross 
J  J 

correlation  of  the  n  seismograms  on  the  interval  of  duration  t 

in  which  e .  is  a  maximum: 

J 


=  e  .  p  peak 

J 


When  G.  is  plotted,  the  shallow  earthquakes  are  emphasized  more 
J 

than  the  events  with  no  pP  in  the  coda  interval.  Typically, 

P  for  shallow  earthquakes  is  .5,  whereas  for  surface  and 

deep  events  it  is  close  to  zero.  (See  Table  1). 

The  test  has  been  applied  to  twelve  earthquakes  and 

six  surface  events.  Plots  of  for  earthquakes  are  shown  in 

J 

Figures  7  and  8.  Plots  for  the  surface  events  are  not  shown, 

as  ii,  for  all  six  was  found  to  be  negligible.  Maximum  values 
J 

of  r|.  for  surface  events  and  earthquakes  are  given  in  Table  1, 

J 

along  with  the  h.*s  at  which  they  occurred.  For  surface  events, 

J 

the  mean  Tlma:X  was  .11  while  for  shallow  earthquakes  it  was  5.1. 

Values  of  depth,  h.,  agree  very  well  with  the  nominal  depths  in 

J 

nine  of  the  ten  cases,  standard  deviation  being  2.8  kilometers. 

The  earthquake  of  15  July  1963,  while  nominally  at  60  km  depth, 
registered  a  distinct  peak  at  19  km  on  the  7  diagram.  (This 
anomaly  will  be  investigated  more  closely.)  The  25  August  1963 
earthquake  is  peculiar  in  that  it's  t)-h  diagram  has  three  distinct 
■peaks,  but  these  may  have  been  caused  by  water-air  and  earth-water 
interfacial  reflections  under  the  Sea  of  Okhotsk.  This  problem 
has  not  been  resolved  conclusively  yet. 

Further  research  will  be  reported  on  this  topic  as 
more  events  are  tested.  Also,  the  use  of  Ppeak  as  a  possible 
discriminant  between  shallow  earthquakes  and  shots  bears  further 
statistical  investigation. 
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APPENDIX  A 


1 .  General  Nonstationary  Case 

Definitions  and  assumptions: 

n(t):  a  zero  mean,  not  necessarily  stationary, 

Gaussian  process,  with  n(t)n( x)  =  R  (t,x) 
x(t):  a  pulse  "signal"  for  which  the  Fourier  Transform, 

X(f)  =  /x(t)e“^2  7rftdt ,  exists 
a(t):  the  "time  window";  it  is  assumed  that  A(f)  exists 

Z ( t )  =  a(t) [x ( t ) +n( t ) ] 

With  these  definitions,  the  following  equations  may  be  written 
by  inspection. 

Z(f)  =  A ( f )  8  X(f)  +  /a(t)n(t)e~  2irft  dt  (A-l) 

Zx(f)  +  Zn(f)  (A-2) 

E[|Z(f)|2]  =  / / dtd ia( t ) a( t) [x(t )+n(t) ] [x(i )+n( t) ]e_j 2 ( t_ t) 

=  |  A  ( f )  &X  ( f )  |  2  +  J/dtdta(t)a(  T)Rn(t,T)e-t527rf(t_T) 

=  | Z  x ( f ) |2  +  |Zn(f) |2  ( A—  3 ) 

2 

The  expected  value  of  |Z(f)|  is  simply  the  sum  of  the  signal 
and  noise  contributions.  It  should  be  noted  for  future  re- 
ference  that  the  integral  defining  | Z^ ( f ) |  must  be  real  and 
nonnegative  for  any  legitimate  R  (t,x). 

covtlZCfj,)  |2,  |Z(f2)  |2]  =  E[|Z(f,)  |  2  |  Z  ( f  2 )  |2]  - 

-  EC | Z(f1 )  | 2 ]E[ | Z ( f 2 )  |2]  ( A—  4 ) 

E[ | Z(f1) |2 |Z(f2) |2]  =  ////dtdTdodea(t)a(T)a(a)a(e)  * 

-^irf^t-x)  -j  2  irf  2  (  a-B ) 

•  e  e 

•  Lx(t)+n(t)][x( i)+n( t) j [x ( a ) +n ( a ) ] [x ( B ) +n ( B ) ]  (A- 5) 
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When  the  integrand  is  multiplied  out  and  expected  values  are 
taken,  only  terms  with  0,  2  or  4  factors  of  n( • )  will  be  non¬ 
zero.  To  simplify  the  expression  involving  all  four  of  the 
n(*)>  the  following  property  of  zero  mean  Gauss  processes  is 
employed . 

n(t)n(x)n(a)n(e)  =  Rn ( t , x ) Rn ( a , 3 ) +Rr ( t , a ) Rr (  x ,  3  )  + 

+  Rn(t,3)Rn( x,a)  ( A-6 ) 

Define 

-j  2  tt(  f  t  +  f  ?a  ) 

F(f1,f2)  =  //dtdaa(t)a(a)Rn(t,a)e  X  (A-7) 

and  note  that  F(f.  ,f2)  =  F(f2,f  )  =  F*(-f  ,-f2)  and  that 
|Zn(f)|2  =  F(f,-f).  With  this  notation,  the  various  contri¬ 
butions  to  the  integral  of  Equation  (A-5)  may  be  tabulated  as 


follows . 

Tex-m  in  Integrand 

Contribution  to  Integral 

1) 

n(t)n(x)n(a)n(3) 

F(f1,-f1)F(f2,-f2)+|F(fi,f2) +|P(f1> 

2) 

x(t)x( x ) x ( a ) x ( 3 ) 

lzx(fl> 1 2 1 Zx  <  r2 ) I2 

3) 

x(t)x( x) n( a ) n ( 3 ) 

|Zx(f1) |2F(f2,-f2) 

4) 

n( t )n( x ) x ( a ) x ( 3) 

|zx(f2)  |2F(f1  ,-q) 

5) 

x( t)x(a)n( x)n( 3) 

zx(fi)zx(r2)F<-f1>-f2) 

6) 

x(t)x( 3 ) n ( x )n( a) 

zx(ri)zx«(f2)F(-fi,f2) 

7) 

x( x)x ( a)n( t )n( 3 ) 

Zx*<fl)Zx<f2)P(fl>-f2) 

8) 

x( x )x( 3 )n( t )n( a) 

Zx*(fl)Zx*(f2)F(f1,f2) 

Since 

1  z  (  t-l  )  |  2  =  | 

Zx(fi)|2  +  F(f1,-f1) .  (A-8) 

E[ | Z ( f 1 ) | ^ ]E [ | Z ( f 2 ) | ^ ]  will  have  four  terms,  which  are  identical 
to  the  terms  appearing  in  lines  2-H  and  the  first  entry  in  line  1 
of  the  above  table.  Thus  the  covariance  (Equation  (A-4))  is  given 
by  the  second  and  third  entries  of  line  1  plus  lines  5-8. 
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cov[ |Z(f1)  |2, |Z(f2)  |2]  =  |F(f1,f2)  |2  +|P(f1,-f2)  | 

+  2Rfe[Zx(f1)[Zx(f2)P(-f1,-f2)+Zx»(f2)P(-f1>f2)]  ( A—  9 ) 

This  specializes  to 

var[|Z(f)|2]  =  |  F  (  f ,  f )  | 2  +  |F(f,-f)|2 

+  2ReZ  2(f)F(-f,-f)+2 |Z(f) |2F(f,-f)  (A-10) 

A.  A 

where  the  fact  that  F(f,-f)  is  real  has  been  used  in  dropping 
the  real  part  operator  from  the  last  term. 

2 .  Stationary  Case,  Rn(t,x)  =  R  (t-x) 

In  the  case  where 

Rn(t,x)  =  Rn(t-x)  =  /Sn(f)eJ27rf(i_  x)df  (A-ll) 

the  expressions  of  the  previous  section  may  be  simplified. 
Specifically , 

-  j  2ir f ,  t  -j2irf?ot 

F ( f x , f 2 )  =  //dtdaa(t)a(a)R  (t-a)e  e 

-j2irf,t  -J2irf_t 

=  /dta(t)e  [a(t)e  £  Rn(t)] 

The  term  in  the  brackets  has  a  Fourier  Transform 

A(f+f2)Sn(f) 

Thus  the  total  integral  is 

F(f1,f2)  =  A ( f )  £  A(f+f2)Sn(f) | f=f 

=  JduA(f1-u)A(u+f2)Sn(u)  =  JduSn(u)A*(u-f1)A(u+f2)  (A-12) 
This  gives 

| Z( f ) | 2  =  |Zx(f)|2  +  JduSn(u) | A(u-f ) | 2  (A-13) 
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var [ | Z ( f ) | 2 ]  =  | /duSn(u) |A(u-f) |2 |  V| JduSn(u)A*(u-f)A(u+f) |2 

+  2 | Z  ( f ) | 2/duS  (u) | A(u-f) +2ReZ  2 ( f ) JduS  (u)A#(u+f)A(u-f ) 

A  II  A  1 1 

(A-14) 

Cov[ | Z ( f^ ) | £ , | Z ( f ^ ) | 2 ]  may  be  obtained  by  substituting  (A-12) 
into  ( A— 9 ) • 

A  special  case  of  some  interest  occurs  when 

A ( f )  =  0  for  | f |  >  W  ( A-15) 

If  (A-15)  holds,  then  for  any  f.,  ,f ,  such  that 


F(f-L,±f2)  will  be  zero  and  thus 

cov[|Z(f1) |2,|Z(f2) |2]  =  0  (A-16) 

Furthermore,  if  (A-15)  holds 

A(u+f)A*(u-f )  =  0  for  | f |  >  W  (A-17) 

and,  therefore,  for  |f|  >  W 

var [ | Z ( f ) | 2 ]  =  |/duSn(u) |A(u-f) |2|2+2|Zx(f) |2/duSn(u) |A(u-f) |2 
=  (  I zn ( f )  I  2 ) 2  +  2  I Zx(f)  I  2 | Zn( f ) |2  (A-18) 

3 .  Nonstationary  Case  Using  n  (t)  =  n(t)-n(t)  as  Derived  in 

Section  F 

In  Section  F,  an  estimate  of  n(t),  n(t),  that  is 

2 

optimal  in  the  sense  of  minimizing  ng  (t),  where 

ne(t)  =  n ( t ) -n ( t )  (A-19) 

and  n(t)  is  based  only  on  values  of  n(t)  for  t<0,  was  derived 
by  considering  the  hypothetical  system  shown  below. 
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w(t ) 


realizable  realizable 

whitener  inverse 

This  estimate  was  given  by 

0 

n(t)  =  /  w ( a ) h ( t-a ) da 

—  oo 

which  gives 

(t) 

n(t)  =  /  w(a)h(t-a)da 
0 

where  the  upper  limit  in  parentheses  is  redundant  since 

is  the  impulse  response  of  a  realizable  system. 

For  convenience,  define 

_ „  (t)  (x)  _____ 

b ( t ,  t  )  =  n  (t )n  (t )  =  /  da/  dgw(a)w(B)h(t-a)h(t- 

0  0 


h  ( t ) 


n(t) 


h  1(t) 


Since 


w( a)w( 6 )  =  Uq ( a-6 ) 


Equation  (A-22)  may  be  simplified  to 


b( t  ,t )  =  { 


I 


(min(t,x ) )dah(t_o)h( T_a) 


Letting 


v  =  min( t , t)  -  a 


this  may  be  rewritten  as 
min( t , t ) 

/  dvh( v)h( v+ | t-x | ) 

b  ( t ,  t  )  =  {  0 


To  contrast  this  with  R  (t-x),  note  that 


Rr(  t  )  =  h(  t)  8  h(  -  t) 


min( t ,  i) >0 
min(t , t) <0 


min( t , t ) >0 
min ( t ,  x)  <0 
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n  ( t ) 

( A-20 ) 

( A-21 ) 
h  ( t ) 

) 

(A-22) 

(A-23) 

( A—  2  ^4 ) 

( A-2  5 ) 

(A-26) 

(A-27) 
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and  therefore 


Rn(t-T) 

Now  defining 


/  dvh(v)h(v+ 1 1  —  t | ) 
0 


(A-28) 


z(t)  =  a(t) [x ( t ) +ng ( t ) ] 

2 

the  moments  of  |Z(f)|  may  be  obtained  by  substituting  b(t,x) 
for  Rn ( 1 3  t )  in  the  expressions  of  Section  1  of  this  appendix. 
This  leads  to 


min(t  ,  t) 

P(fl,f2)  =  //dtdxa( t )a( t) /  dvh( v)h( v+ | t- t | ) • 

—  j  2lT  (  f  ,  t  +  f  „  T  ) 

•  e  1  *  (A-29) 


Separating  this  integral  into  the  two  regions,  t < t  and  t>x 
and  substituting 

a  =  min(t,x),  B  =  1 1 —  x [  (A-30) 

leads  to 

°°  t  a  -j  2 XT  [f,  (  a+B  )+f  ?a] 

F(f13f2)  =  /  dt /  dx a( a) a( a+B ) /  dvh ( v ) h ( v+B ) e 

—  00  -  00  Q 

00  “  a  -  j  2  xr  [  f ,  a+f  (  a+  B  )  ] 

+  /  dt /  dxa( a)a( a+B) /  dvh ( v) h( v+B ) e 


-oo  t 

0 

(A-31) 

but 

»  t 

/  dt/  dx  = 

00  oo  00  00 

/  dx /  dt  =  /  da/  dB 

-oo  X  —  oo  0 

OO  00 

(A-32) 

and 

oo  oo 

/  dt /  dx  = 

/  da/  dB 

( A—  33) 

_oo  t  -00  0 


Therefore , 


00  00 


F( f ,  j  f p )  =  /  da/  d $J  dva( a) a( a+B )h( v)h ( v+B ) e 


-j2n(f1+f2)o 


0  0 


-  j  2 TT f  1  B  -j2Trf2B 

(e  +  e  ) 


( A- 3  4 ) 


and,  in  particular, 

00  00  £ 

F  (  f ,  -  f )  =  2  /  da/  dB  /  dva(  a )  a(  a+B  )  h  (  v)  h  (  v+B  )  cos2  irf  B  (A-35) 

0  0 
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F(f,+f) 


2 /  da/  d6 /  dva(a)a(a+6)h(v)h(v+B)e  ^ 

-°°00  ( A—  36 ) 


To  convert  Equations  (A-3^)  to  (A-36)  to  the  stationary 
case  (for  which  the  autocorrelation  function  is  given  by  Equa¬ 
tion  (A-28)  rather  than  (A-26))  it  is  only  necessary  to  change 
the  upper  limit  on  the  integration  from  a  to  +<*>.  With  this 
change.  Equations  ( A-3*0 -( A-36 )  can  be  simplified  as  has  already 
been  demonstrated  by  Equation  (A-12).  Without  this  change,  it 
does  not  appear  possible  to  simplify  these  equations  further 
without  making  specific  assumptions  about  a(t)  and  h(t).  One 
set  of  assumptions,  which  leads  to  some  simplification, 
is  considered  in  a  later  section  of  this  appendix. 


4.  Choice  of  n(t)  to  Minimize  |Z  „(f) 

ne 


Considering  now  an  estimate,  n(t),  of  the  form 
0 

n(t)  =  /  g(t , a ) n ( a ) da 


( A—  37 ) 


the  problem  is  to  choose  g(t,a)  to  minimize  |Zne(f)|  .  [Note, 
by  inspection  of  Equation  (63),  it  can  be  seen  that  the  choice 
of  n(t)  considered  in  the  previous  section  corresponds  to 


g  ( t ,  a )  =  /  duh(t-u)h  ‘'‘(u-a)  ] 

(  a) 


Substituting 


0 

n  (t)  =  n(t)  -  /  g(t,a)n(a)da 


(A-38) 


(A-39) 


into  the  integral  defining  |Z  (f)|  and  taking  expectations 

ne 

yields 

0 


zn  (f)|2  =  //dtdia(t)a(  x)e  J’2lTf(t  x){Rn(t-x)-/  deg(x,e) 

6  —  00 

0  0 

- R  < t  —  6 ) — /  dag(t,a)R  ( x-a )  +  / /dadgg(t,a)g( t  ,  g ) 


Rn(a-e) } 


(A-40) 
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Letting 

g ( t , a )  =  gQ ( t , a )  +  eg£ ( t ,a)  (A-4l) 

the  problem  is  to  choose  gn(t,a)  in  such  a  way  that  | Zn  (f)| 

u  e 

is  a  minimum  at  e=0  for  any  g  (t,a). 


ll|Z 


(f)|2  =  //dtdia(t)a(  i)e  :  2  "f  +  t){0 -/ 


e  =  0 


dBg£(T ,8) 


0  0 

•Rn(t-6)-/  dag£(t Ja)Rn(T-a)+{/dad6[g0(t,a)gt ( t,8)  + 


+S0(T , B)ge(t , a) ]Rn( a- 6)  } 

=  //dtdTa(t)a(  T)e-J<::  ^  {/  d8ge(Tj8)[/  dagQ(t,a)- 

—  00  —  oo 

0  0 

■Rn(a-8)-Rn(t-B)]  +  /  dag£(t ,a) [/  d8gQ( t , 6 )Rn( 8-a)  - 

—  OO  —  oo 

-Rn(T-a)]}  (A-42) 

A  sufficient  condition  for  Equation  (A-42)  to  be  zero,  for  any 
choice  of  g  (t,a)  is 
0 

/  dagQ(t 5a)Rn(a-6)  =  R  (t-8)  for  8<0  (A-^3) 

—  OO 

Guessing  that  the  solution  to  the  integral  equation  is  given  by 
Equation  (A-38),  which  corresponds  to  the  earlier  derivation 
of  n( t ) ,  yields 

gQ ( t , a )  =  / duh( t-u)h-1(u-a)  (A-44) 

/  dagQ(t  ,oc)Rn(  a-6)  =  /  da/duh ( t-u ) h  1 ( u- a ) Rn ( a-8 ) 

—  OO  —00 

=  duh  ( t-u ) /doth  ^(u-a)R  (a-8) 

—  oo 

0  -1 

=  /  duh( t-u) [h  ^(u)  8  Rn( u- 8 ) ] 

—  00 

[equation  continued] 
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0 

=  /  duh( t-u)h[-( u-8 ) ] 

—  cr 

Rn(t-B)  for  8<0 

^  min( 1 ,  8 ) 

R  ( t —  6 ) — /  h(t-u)h( 8-u)du  for  8>0 

0 

=  Rn(t-B)-b(t,6)  (A-45) 

where,  of  course,  b(t,8)  is  zero  for  8<0.  Thus  g|(t,a),  as 
given  by  Equation  (A- 4  4)  is  a  solution  to  Equation  (  A—  ^4  3 )  - 
To  show  that  this  does  indeed  give  a  minimum  of 
|zn  (f)|  ,  it  is  necessary  to  show  that  the  second  derivative  is 
positive  at  e  =0. 

— p|Z  (f)|^  =  f  fdtdra(t)  a(  T)e~j  ^7rr(t~T  '  f  da/  dgg.(t,a) 

9e  e  -oo  -oo 

*g0 ( x, 8 ) Rn( a- 6 ) 

=  //dtdxa(t)a(  i)e-t  2  T'>  /  dBgQ  (  t  ,  8  )  Rr(  t- 6  ) 

—  00 

=  //dtdia(t)a(  i)e"J2lTf(T'_T)[Rn(T-t)-b(  T,  t)]  (A-46) 

where  Equation  ( A— 4 3 )  is  used  for  the  first  integration  and 
Equation  (A-4S)  for  the  second. 

Since  n(t)  and  n  (t)  are  independent,  zero  mean  processes 

and 

n  ( t )  =  n(t)  +  ng(t)  ( A—  ■47) 

it  follows  that 

Rn(t-x)  =  n(t)n(t)  +  b(t,i  )  (A-48) 

and  therefore  the  quantity  in  the  brackets  in  Equation  (A-46)  is 
the  autocorrelation  function  of  n(t).  Therefore  Equation  (A-46) 
is  simply  |Z^(f)|  ,  which  cannot  be  negative. 
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5 .  Nonstationary  Case  with  Rectangular  a(t)  and  h(t) 

In  an  effort  to  obtain  some  measure  of  the  improvement 
attainable  by  using  n(t),  specific  assumptions  about  a(t)  and 
h(t)  are  now  considered.  Unfortunatley ,  even  with  these  assump¬ 
tions,  only  rather  limited  results  have  been  obtained. 

Specifically,  it  is  assumed  that 


1  0  <t  <T 

a(t)  =  { 

0  elsewhere 

1  0<t<B 

h  ( t )  =  { 

0  elsewhere 

From  this  it  is  easy  to  show 


(A-49) 


(A-50) 


A(f)  =  e"J27rf(T/2)  sln*fT 


(A-51) 


Sn(f)  =  | H( f ) 


,  sin  irfB  >  2 

7T  f 


(A-52) 


Pursuing  these  assumptions, 

00  oo 

F(f.-f)  =  2/  dB  /  daa(a)a(a+B) 


T-B  for  B<T 


a 

/  dvh(  v)h(  v+B  )  cos2  Ttf  B 


min( a ,B-B ) 


If  B>T 


T  T-B  T  , 

2  j  dBcos27i  f  B  /  da*a  =  /  dBcos27if  B  (  T-B  ) 

0  0  0 

l/SimrfTN2  ^  /sin7ifTN2 
2( - M  (  rf  } 


( A-53 ) 


If  B<T 

B  B-B  T-B 

F(  fj  -f )  =  2/  dBcos2irfB[/  da-a  +  /  da(B-B)] 
0  0  B-B 


[equation  continued] 
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=  /  d3cos27rfg[(B-3)"  +  2(B-3)  (T-B)  ] 

0 

_  l,sinirfB,2  ,  sinirfBx  2  , /m  ^/SimrfB^ 

-  2"( - —)  «  ( - —)  +(T-B)( - itf- ) 


(A-54) 


In  summary 


Z  <f) 
e 


=F( f ,  -f )  =  { 


1 , sinufT, 2 

2  Ttf 

1/  sin  7rfBN  2, 
2"  ttT  * 


/Sin  7rfT>.  2 
1  77  f 


B>T 


(A-55) 

( -1— "  )  2  +  (  T-B )  ( — b )  2  B<T 

7TI  ~  + 


irf 


These  should  be  compared  to  the  corresponding  expression  for  the 
same  a(t)  and  h(t),  but  without  n(t).  This  is 


Z  (f) 
n 


-  Sn<f> 


a  I A  ( f ) 


,  simrfTv  2 
'  Ttf  ' 


^  simrfB ^  2 
7T  f 


(A-56) 


Examining  a  few  special  cases. 


i)  for  B=T  \Zn  (f) |2 
n< 


l\r,  /(,N  |2  1  /  sinirfT  N  2  ,sinirfTN2 

2'Zn(fJl  =  2(  *f  }  8  (  fff  } 

(A-57) 


ii)  for  B> >T  B(sin!fT)2  «  | ( f )  | 2  >>  |Z  ( f )  | 2 

TTi  1  n  1  1  n  * 

e 

_  1  /  simrfT  N  2  „  /SinwfTN2  ,  „  r0N 

-  2(  '  Ttf  }  8  (”  Trf  }  (A"58) 

iii)  for  B«T  |  ZR(  f )  |  2  -  |Z  ( f )  |  2  *  TC-1^-8)2  (A-59) 

e 

The  equations  indicate  an  "improvement"  of  1/2  when 
the  maximum  correlation  time  of  the  noise  equals  the  length  of 
the  observation  interval,  a  large  improvement  when  the  interval 
is  much  shorter  than  the  maximum  correlation  time  of  the  noise, 
and  essentially  no  improvement  if  the  correlation  time  of  the 
noise  is  short  compared  to  the  total  observation  interval. 
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APPENDIX  B 


Weighting  Factors 

Assume  that  s^(t)  +  n^(t)  and  SgCt)  +  ng(t)  are  two 
seismograms  of  different  signal-to-noise  ratio  which  are  to  be 
added  together  coherently.  A  weight  or  bias  a  is  to  be  given 
to  the  second  so  that  the  signal-to-noise  ratio  of  the  sum 
seismogram  is  a  maximum: 


=  Maximum 


It  is  appropriate  in  this  case  to  assume  that  the  noise-cross- 
noise  term 


as  stations  1  and  2  are  chosen  sufficiently  far  apart  for  two 
noises  to  be  independent,  and  the  noise  is  assumed  to  be  Gaussian. 


Set 


Then  our  expression  becomes : 


*  All  integrals  are  assumed  to  be  over  unit  time  for  expressional 
clarity. 
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d.  2  2  r 

S1  +  s2  a  +  2aJ  si(t)s2(t)dt 

2  2  2 
+  a  n2 


Since  the  signals  s-^(t)  and  s2(t)  are  coherent,  i.e, 
mutually  proportional,  the  integral  may  be  written: 


J  o1(t)s2(t)dt  =^Js12(t)dt  J  s22(t)dt)‘; 


-  b1b2 

(because  the  Schwartz  inequality  reduces  to  an  equality  in  this 
case.)  The  ratio  to  be  maximized  then  becomes 

(s1  +  as2)2 

2  T  2 

n^  +  a  n2 

For  ease  in  differentiation,  we  make  the  substitution 


(si/s2 )  -  k-j_  (ni^n2)  ~  k2 


( sx  +  as2) 

- 2 - T*" 

n^  +  a  n2 


(k^  +  a) 

2  2 
k2  +  a41 


2 


s. 


n. 


To  find  the  maximum  with  respect  to  a,  set 

2  2 


d_  (■(kl+a>  \  ‘2  .o 
5a  v  k2+a2  J 


or 


2  2 

k9  +  a  -a(k,  4a)  =  0 

^  2  1 

k2  -ak-j^  =  0 
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9 


a 


The  sum  seismogram  is  therefore 


2 

ni 

S(t)  =  s1(t)  +  r1(t)  +  -j- 


2  / 

— 2  (^32^t^+n2^tV 
n2 


Multiplying  S(t)  +  N(t)  by  a  constant  does  not  change  the  signal- 
to-noise  ratio,  so  we  may  write  an  equivalent  sum  as 


f 


So(t) 


Generalizing  this  to  an  arbitrary  number  of  seismograms,  we  find 
that  the  proper  weighting  factor  (for  seismogram  i)  to  maximize 
total  signal-to-total-noise  ratio  is 


ni  jni2(t)dt 


Note : 


In  practice,  signal-to-noise  ratio  is  much  greater  than 

one,  so  that  we  can  approximate  the  true  signal  power 

2 

Jsi2(t)dt  by  J(si(t)  +  n± ( t ) )  dt. 
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b)  the  theoretical  aspects  of  clustering  seismometers  in  a 
large  aperture  seismic  array, 

c)  the  results  of  tests  of  a  revised  automatic  pP  phase 
detection  technique. 

A  method  has  been  devised  to  estimate  the  energy  density  spectra 
of  transient  events  in  such  a  way  that  the  variance  of  the  spectral 
estimates  produced  by  additive  noise  is  minimized. 

The  study  of  clusters  shows: 

1)  dependence  of  array  S/N  performance  on  cluster  spacing, 

2)  a  way  in  which  to  determine  the  amount  of  clustering  possi¬ 
ble  without  S/N  degradation  depending  on  noise  conditions 
at  the  receiving  site. 

Tests  with  the  revised  automatic  pP  detector  have  shown  100$ 
correct  results  with  surface  focus  events  (6),  83$  correct  results 
for  earthquakes  (12),  with  the  remaining  earthquakes  producing 
anomalies . 
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